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Ph Abstract 



In this article, we describe an approach for solving partial differential equations with 

general boundary conditions imposed on arbitrarily shaped boundaries. A contin- 

I— I uous function, the domain parameter, is used to modify the original differential 

equations such that the equations are solved in the region where a domain param- 

^ eter takes a specified value while boundary conditions are imposed on the region 

^~~^ where the value of the domain parameter varies smoothly across a short distance. 

S^ The mathematical derivations are straightforward and generically applicable to a 

If^ wide variety of partial differential equations. To demonstrate the general applicabil- 

r~^ ity of the approach, we provide four examples herein: (1) the diffusion equation with 

<^ both Neumann and Dirichlet boundary conditions; (2) the diffusion equation with 

both surface diffusion and reaction; (3) the mechanical equilibrium equation; and (4) 

• • the equation for phase transformation with the presence of additional boundaries. 

. ^H The solutions for several of these cases are validated against corresponding analyt- 

/\ ical and semi-analytical solutions. The potential of the approach is demonstrated 

3 with five applications: surface-reaction-diffusion kinetics with a complex geometry, 

Kirkendall-effect-induced deformation, thermal stress in a complex geometry, phase 

transformations affected by substrate surfaces, and a self-propelled droplet. 
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1 Introduction 



The smoothed boundary method pl2ll3] has recently been demonstrated as 
a powerful tool for solving diffusion equations with no-flux boundary condi- 
tions imposed at irregular boundaries within the computational domain. The 
method's origin can be traced to the embedded boundary method and the 
immersed interface (boundary) method, both of which embed a more compli- 
cated domain in a computational box with simpler geometry. These methods 
are advantageous because they eliminate the need for a structural mesh when 
solving partial differential equations within the embedded geometries because 
the grid system is obtained by a discretization of the regular computational 
box. (For an overview, see Refs. JUSIEIEIHEJ-) To impose boundary conditions 
at the immersed interfaces, a discretized Dirac delta function is employed 
to distribute a singular source over nearby grid points. Various studies have 
examined optimal discretization of the Dirac delta function [7|10llll] . Simi- 
larly, the level set method can also be considered an immersed-interface-type 
method because the boundary defined by the contour of the zero level set is 
embedded within a regular computational box. Although the level set method 
was developed mainly for tracking moving boundaries [12J , it is also applicable 
for solving partial differential equations with boundary conditions imposed at 
the zero-level-set contour using a technique similar to the immersed interface 
method P^|IT^ . In addition to the methods above, the phase field approach 
possesses certain similarities to embedding interfaces within the computational 
box and also has the significant advantage of avoiding the need to explicitly 
track the interfaces. However, phase field methods are not widely employed in 
simulations that involve explicit boundary conditions along interfaces. While 
the Gibbs-Thomson boundary condition is automatically imposed in the stan- 
dard phase field model, there have been only a few studies in which phase 
field models were used with explicit boundary conditions at interfaces. For 
example, in solidification problems, equilibrium conditions, such as equilib- 
rium temperature or concentration [T5p6] . are imposed at solid-liquid inter- 
faces in which the order parameter field and the temperature field are coupled 
via a latent heat term. Except for this type of phase field model, the direct 
application of boundary conditions at interfaces is rarely used because the 
construction of boundary conditions requires the tedious process of formally 
including an additional energy term in the energy functional, as suggested by 
Cahn [17J. Examples of imposed boundary conditions in the phase field model 
using modified energy functionals can be found in the recent works of Warren 

et ai. [rgfig] . 

In contrast to the techniques for distributing a singular source of boundary 
conditions to grid points near the interfaces in the immersed interface method, 
the smoothed boundary method spreads the zero-thickness boundary into a 
finite-thickness diffuse interface using a phase-field- like, continuously transi- 



tioning domain indicator function (hereinafter termed "the domain parame- 
ter"). Mathematically, this method approximates a Heaviside step function 
as a hyperbolic tangent function having one specified uniform value in a do- 
main and continuously changing its value across the interface to another value 
specifying the other domain. Therefore, boundary conditions are straightfor- 
wardly distributed among the grid points residing within the interfacial re- 
gions in which the domain parameter varies smoothly across a short distance. 
This method has been successfully employed in simulating diffusion processes 
[20pT] and wave propagation [Tp|3l)22f23] constrained within geometries de- 
scribed by a domain parameter with a no-flux boundary condition imposed 
on the diffuse interfaces. Similar approaches have also been proposed to solve 
differential equations constrained in domains defined by order parameters in 
the phase field model P^25f26] . These works demonstrated the potential for 
this type of numerical method that circumvents the difficulties associated with 
constructing a finite element mesh (e.g., meshing the surface and then build- 
ing a volumetric mesh based on the surface mesh or combining regular sub- 
domains that can be easily meshed). Such an approach is particularly useful 
when complex structures are involved. However, the method was only applica- 
ble to no-fiux boundary conditions, and no further extensions to other types 
of equations or boundary conditions have been reported. Recently, Lowengrub 
and coworkers |27ll28p29f30ll3ip32f33j developed an alternative formulation for 
solving partial differential equations with various boundary conditions, based 
on asymptotic analyses commonly conducted in phase field modeling, which is 
different from the general derivation of the smoothed boundary method pre- 
sented in this paper. Although such an implementation for imposing boundary 
conditions differs from the 'formal' practice suggested by Cahn [TTj, it dra- 
matically simplifies the formulation, provides a justification of the method, 
and increases the applicability of the approach. 



In this study, we provide a mathematically consistent smoothed boundary 
method and a precise derivation for the equations, such that the method is 
generalized from its limited original application to a wide range of differential 
equations and boundary conditions. We consider the following specific equa- 
tions: (1) the diffusion equation with Neumann and/or Dirichlet boundary 
conditions; (2) the bulk diffusion equation coupled with surface diffusion and 
reaction; (3) the mechanical equilibrium equation for linear elasticity; and (4) 
the Allen-Cahn or Cahn-Hilliard equations with contact angles as boundary 
conditions. The method is especially useful for three-dimensional image-based 
simulations because of its efficiency and flexibility in handling complex geome- 
tries without structural-mesh techniques. 



2 Background 



The method is based on a diffuse interface description of different phases, 
similar to the continuously transitioning order parameters in the phase field 
method |34|l35|36|37|l38|39j often employed in simulating phase transforma- 
tions and microstructural evolutions in materials. Phase field models are based 
on thermodynamics and kinetics of multiphase system, in which phases (e.g., 
hquid, solid, vapor, or two solids or liquids with different compositions) are de- 
scribed by one or more order parameters with prescribed bulk values for each 
phase. At the interface, the order parameter changes in a controlled manner. 
Asymptotic analyses [HH] can be used to show that the phase-field governing 
equations approach the corresponding sharp interface equations in the sharp 
interface limit. 



Despite the advantages of phase-field-type diffuse interface methods for front 
tracking problems, we focus on another important advantage for efficiently 
solving differential equations within diffuse-interface-defined domains. Here, 
we adopt the concept to describe internal domain boundaries with an order- 
parameter-like domain parameter, which may or may not be stationary and 
takes a value of 1 inside the domain of interest and outside. The equa- 
tions are solved where the domain parameter is 1, with boundary conditions 
imposed where the domain parameter is at an intermediate value (approxi- 
mately 0.5). Figure [I] schematically illustrates the sharp and diffuse interfaces. 
In the conventional sharp interface description, the domain of interest is Q and 
is bounded by a zero-thickness boundary, denoted dD.; see Fig. [lla). Within fi, 
the partial differential equations are solved according to the boundary condi- 
tions imposed at dfl. Conversely, in the diffuse interface description, we employ 
a continuous domain parameter, which is uniformly 1 within the domain of 
interest and uniformly outside. In this case, the originally sharp domain 
boundary is smoothed to yield a diffuse interface with a finite thickness given 
hj < ip < 1. The system thus determines the boundary by variation of the 
domain parameter. In addition, the gradient of the domain parameter V^ 
automatically determines the inward normal vector of the contour level sets 
of ip; see Fig. lite). Our goal is to solve partial differential equations within 
the region where ip = 1 while imposing boundary conditions at the narrow 
transitioning interfacial region where < ip < 1. However, the convention 
can be reversed such that the domain is defined hj ip = 0, in which case the 
following derivation could be modified by replacing ip with 1 — ip accordingly. 
This could be used to solve a problem where multiple equations govern differ- 
ent regions within the computational domain. Furthermore, these equations 
can be coupled through the shared boundary conditions, making the method 
highly versatile. 



3 Formulation 



3.1 General Approach 



The general approach is as follows. The domain parameter describes the do- 
main of interest {ip = 1 inside the domain, and ip = Q outside). The transition 
between the two values described is must be smooth so that the gradient is well 
defined. In this work, we have assumed the domain parameter to take the form 
of a hyperbolic tangent function for three reasons. First, it can be numerically 
implemented with ease. Second, it is consistent with the solution to the phase 
field equations, and thus this choice allows coupling of the two approaches; 
that is, one could for example simulate microstructural evolution using the 
phase field model, but impose interfacial flux boundary condition using the 
smoothed boundary method. Third, when given a non-smooth microstructural 
data, one could use the phase field equations to obtain the domain parameter 
from the discontinuous data. Other forms of domain parameters are possible, 
as long as the transition is monotonic and the gradient of the domain param- 
eter has a narrow peak in the interfacial region. A better convergence could 
possibly be obtained using a function that has a more confined gradient; how- 
ever, examining other forms of the domain parameter is beyond the scope of 
this paper. As an example, we consider the Laplacian of the function, H. As 
the first step in deriving formulation for the Neumann boundary condition, 
V^H is multiplied by the domain parameter, ip. Using identities of the product 
rule of differentiation such as: 

i,V^H = V ■ i^VH) - V^ ■ VH, (1) 

we obtain terms proportional to Vip. Because the inward unit normal of the 
boundary (pointing to the regions where ip = 1), n, is given by Vip/lVipl, such 
terms can be written in terms of dH/dn = VH -n = VH -Vip/lVipl, and thus 
reformulated to be the Neumann boundary condition imposed on the diffuse 
interface. 

Similarly, to derive the smoothed boundary formulation for the Dirichlet bound- 
ary condition, the differential equation is multiplied by the square of the 
domain parameter. Again using mathematical identities, ip'^V'^H = ipV ■ 

{ipVH) - ipVip ■ VH, where ipVip ■ VH = Vip-V {ipH) - H \Vipf, we obtain: 

iP^V^H = tpv ■ itpVH) -[Vtp-V itpH) - H\V^PW (2) 

Note that H , associated with | Vt/'p appearing in the last term, is the boundary 
value, H\qq, imposed on the diffuse interface. Specific details of the derivation 
depend on the equation to which the approach is applied, and we therefore 
provide four examples below. 



3.2 Diffusion Equation 



The first example is the diffusion equation with Neumann and/or Dirichlet 
boundary conditions. The Neumann boundary condition is appropriate, for 
example, as the no-fiux boundary condition, whereas the Dirichlet boundary 
condition is necessary when the diffusion equation is solved with a fixed con- 
centration on the boundaries. For Pick's Second Law of diffusion, the original 
governing equation is expressed as: 

dC 

— = - V ■ / + ^ = V ■ (DVC) + S, (3) 

where j is the flux vector, D is the diffusion coefficient, C is the concentration, 
S is the source term, and t is time. Instead of directly solving the diffusion 
equation, we multiply Eq. ^ by ijj, the domain parameter that describes 
the domain in which diffusion occurs, and use the identity ipV ■ (DVC) = 
V ■ (ipDVC) — VV' ■ (-DVC) to obtain the smoothed boundary formulated 
diffusion equation: 

dC 

ij— = V ■ i^DVC) - Vtp ■ {DVC) + i)S. (4) 

Next, we consider the boundary condition in this formulation. The Neumann 
boundary condition is the inward flux across the domain boundary, mathe- 
matically the normal gradient of C at the diffuse interface, and is treated 
as: 

where n = V'^p/\V^p\ is the unit inward normal vector at the boundary defined 
in the diffuse interface description. Note that the flux at the interface is equal 
to DBm- Equation (|5| is rearranged to become V?/' ■ {DVC) = — \Vip\ DBjy 
and substituted back into Eq. Q; thus, we obtain: 

^ = iv ■ {i^DVC) + ^^DB, + S, (6) 

with the Neumann boundary condition appearing in the second term. When 
a no-flux boundary is imposed, the second term vanishes and the resulting 
equation is the same as that proposed in Refs. P]I2|3|20II22] . 

To impose the Dirichlet boundary condition, we manipulate the original gov- 
erning equation in a procedure similar to the derivation of Eq. ([6]). Multiplying 
Eq. Q by ^ and using the identity ipVip ■ {DVC) =D[V'^-V {i/jC) - CVijJ ■ 
V^ip] = DlVip ■ V {ipC) — C\Vip\^] to replace the second term, we obtain: 

-V ■ {ijDVC) - ^D[V^Ij ■ V (^C) - Bd |V^|'] + S, (7) 



dt ijj ijj 



where B^ is the Dirichlet boundary condition imposed at the diffuse interface 
to replace C, associated with \Vip\ in the third term. The convergence to the 
imposed Neumann and Dirichlet boundary conditions is shown in Appendix 

El 

In this method, the boundary gradient, -Bat, and the boundary value, B^, are 
not specified to be constant values; they can vary spatially and/or temporally 
or be functions of C or other parameters. In addition, it is convenient to use 
weighting factors to combine Eqs. (Ig]) and ([T]) to impose Neumann and Dirich- 
let boundary conditions simultaneously to yield mixed (or Robin) boundary 
conditions. The equation then becomes: 

^ = iv ■ i^DVC) + ^-^DBnWn ~^W- ^ii^C) - BdW\^]Wd + S, 

(8) 
where W^ and Wd are the spatially dependent weighting factors for the 
Neumann and Dirichlet boundary conditions, respectively {Wn + Wd = 1). 
These factors can be a linear combination when imposing Robin boundary 
conditions or be employed to impose Neumann and Dirichlet boundary con- 
ditions at different regions of the interface. Moreover, a small nonzero value 
(10~^^ < V < 10~^) should be added to the domain parameter appearing in 
the denominators to avoid singularities resulting from the terms 1/ip and l/ip"^ 
in regions where ^ = 0. 



3.3 Coupled Surface-Bulk Diffusion 



The second example demonstrates that surface diffusion can be incorporated 
into the smoothed boundary equation derived above. For this case, we take 
the set of equations that includes the surface reaction, surface diffusion and 
bulk diffusion to describe an oxygen reduction model in a solid oxide fuel cell 
(SOFC) cathode [IQ]. The oxygen- vacancy concentration, C, on the cathode 
surface is governed by Pick's Second Law for surface diffusion: 

where n is the coordinate along the inward unit normal vector of the surface 
and Vg is the surface Laplacian. The parameters Dy,, k, D^, and L are the 
bulk diffusivity, reaction rate, surface diffusivity, and accumulation coefficient, 
respectively. Thus, the term on the left-hand side represents the flux from the 
bulk, and the terms on the right-hand side represent the surface reaction, 
surface diffusion and surface concentration accumulation, respectively |10]. 
Here, these parameters are all assumed to be constant for simplicity. In the 
bulk of the cathode material, the oxygen-vacancy diffusion is governed by 



Pick's Second Law for bulk diffusion: 



dC 

In 



D^V^C. 



(10) 



To simulate the oxygen-vacancy concentration evolution in the cathode, the 
two diffusion equations, Eqs. ^ and (10), are coupled and must be solved 
simultaneously, with the flux normal to the cathode surface as a common 
boundary condition. Recently, a similar set of equations was formulated using 
anothor diffuse interface approach combined with an asymptotic analysis [30] , 
which led to two differential equations coupled by a common boundary condi- 
tion. We show below that we can eliminate the need for solving two separate 
equations by applying the smoothed boundary formulation described herein 
to obtain a single equation that governs both surface and bulk effects. 



Similar to the derivation of Eq. (|6]), we first multiply Eq. (10) by if) and 
apply the product rule of differentiation to obtain the bulk diffusion equation 
containing a boundary term, DbV^ ■ VC, similar to Eq. (|4]). As in Eq. (|5|, 
the normal gradient at the diffuse interface is defined by dC jdn = — VC ■ 
V'0/|V'?/'|. Substituting this relationship back into Eq. ^ and rearranging 
terms gives: 



Vz/^- VC 






kC - D^V^C + L 



dC 

~dt 



(11) 



Using this relation, we replace the boundary term to obtain the smoothed 
boundary formulated equation: 



dC 
~dt 



DbV ■ (V^VC) 






kC - D^VlC + L 



dC 
~dt 



(12) 



which combines the bulk diffusion and surface diffusion terms into a single 
equation used in the examples presented in Sections 4.2 and 5.1 Again, a 



small nonzero value should be added to the domain parameter appearing in 
the denominators. In the bulk (jV'?/'! = and ip = 1), Eq. (12) reduces to 
Eq. (10). When the interfacial thickness approaches zero, it reduces to Eq. ([9]) 
at the interface (|V^| 7^ 0), as proven in Appendix A 



The surface Laplacian (V^ 
gradient given by: 



Vs ■ Vs) is calculated according to the surface 



(I 



n 



n 



Vip Vip 



|v^| iv^L 



(13) 



where I is the unity tensor, '®' is the Dyadic product, and n is the inward 
unit normal vector of the diffuse interface, as used in Ref. [30]. (In indicial 
notation, the surface gradient is expressed as: {5ij — ninj)d/dxj, where 6^ is 



Jij 



the Kronecker delta {6ij = 1 Hi = j, and 6ij = if i 7^ j). The repeated indices 
indicate summation over the index. See Appendix [B] for details.) To simulate 



only surface diffusion on a diffuse-interface-described geometry, one can simply 
eliminate the bulk-related and reaction terms in Eq. ^ to obtain L{dC/dt) = 
-DsVgC, such that the concentration evolves only over the interfacial region. 



3.4 Mechanical Equilibrium 



The smoothed boundary method can also be applied to the mechanical equilib- 
rium equation. When a solid body is in mechanical equilibrium, all forces act- 
ing on the body are balanced in all directions, as represented by daij/dxj = 0, 
where 0"^ is a stress tensor component, which is the force per unit area along 
jth axis on the surface whose normal vector is along the ith axis. Repeated in- 
dices indicate summation over the index. For linear elasticity, the stress tensor 
is given by the generalized form of Hooke's Law: aij = Cijki^Ski — pSki), where 
Cijki is the elastic constant tensor and p is a scalar body force, such as ther- 
mal expansion (aAT) or misfit eigenstrain ((a^ — am)/cim, where Op and am 
are the lattice constants of the precipitate and matrix phases, respectively), 
depending on the governing physics. The total strain tensor is defined by the 
gradients of displacements as Eij = [{dui/dxj) + {duj/dxi)]/2, where Ui is the 
infinitesimal displacement in the ith direction. Substituting Hooke's Law and 
the total strain back into the mechanical equilibrium equation gives: 



d 1 fdtik duA _ _d_ 

dxj ^^ 2 \dxi dxkj dxj 



pCijkiSki ■ 



(14) 



Multiplying Eq. (14) by the domain parameter that distinguishes the elastic 



solid region (ifj = 1) from the environment {^p = 0) and using the product rule 
of differentiation yields the smoothed boundary formulation: 




P5^ 



kl 



_d_ 

dxj 



■ippCijkiSki 



(15) 



see Appendix O for details of the derivation. 



The traction exerted on the solid surface is defined by Ni = —aijUj, where 



n 



Wip/\W^jj\ is the inward unit normal of the solid surface. (In indicial 
notation, dip/dxi = Vip and J{dip/dxi){dip/dxi) = \Vil)\.) Therefore, the 
surface traction force is given by: 



N, = -{C, 



ijkl 



1 fduk dui " 

2 \dxi dxk_ 



P5. 



kl 




(16) 



Substituting Eq. (16) into Eq. (15) yields the smoothed boundary formulated 



mechanical equilibrium equation with a traction boundary condition on the 
solid surface: 



dxj 



^C.,H^ \ dx 



duk dui 



dxk. 



\vm^ 



dXj 



'ippCijkiSki 



(17) 



where d{ip pCijki^ki) / dxj = pt can be treated as an effective body force in the 
ith direction. 



For linear elasticity problems with prescribed displacements at the solid sur- 
face, one can perform the smoothed boundary formulation, as in the derivation 



of the Dirichlet boundary condition in Section 3.2, by multiplying Eq. (14) by 
ip"^ and using the product rule to obtain: 



^ 



dx 



3 L 



"ipCijki 




dip dip 

Uk^ \-ui- — 

OXi OXk 




f dipUk (hpvA 



2\ dxi 



dxk J 



where the displacements Uk and ui appearing in the third term on the left- 
hand side should be the boundary values of the displacements at the solid 
surface; see Appendix [C] for the derivation. A similar formulation for solving 
the mechanical equilibrium equation within a domain defined by a phase- 
field-like order parameter can also be obtained by the asymptotic approach 
by matching terms of different orders j41j . 



3. 5 Equations for Phase Transformations with Additional Boundaries 



Phase transformations affected by a mobile or immobile surface or other 
boundaries are of importance in many materials processes, including hetero- 
geneous nucleation that occurs at material interfaces P^JITU] . Maintaining a 
proper contact angle at the three-phase boundary (where the interface be- 
tween the two phases meets the surface) is necessary to capture the dynamics 
accurately because the contact angle represents the difference between the sur- 
face energies (tensions) of the different phase boundaries. Although researchers 
have previously developed methods for imposing contact-angle boundary con- 
ditions on sharp domain walls [IHIIS], here we show that a similar model 
with diffuse domain walls can be obtained simply by applying the approach 
described above. Below, we assume that the boundary is immobile, but this 
assumption can be easily removed by describing the evolution of the domain 
parameter as dictated by the physics of the system. 

In the AUen-Cahn and Cahn-Hilliard equations of the phase field model, the 



10 



total free energy has the following form [M1I35] : 

/(0) + YI^0I' dn, (19) 

where is the phase field order parameter commonly used to define differ- 
ent phases, /(</>) is a double-well free energy functional (in terms of 0), e 
is the gradient energy coefficient, and Q is the domain of interest. At the 
extremum of the functional F, the variational derivative of the total free en- 
ergy vanishes: 6F/6(j) = 0. This requirement provides the following conditions: 
df/d(j) - e^VV = G fi, which can be reformulated as V/ = W{e^\W(p\^)/2, 
by multiplying both sides by V0. We thus find a useful equahty for deriving 
the contact angle boundary condition: |V0| = y/2f/e; see Appendix |d] for 
details. 

In the smoothed boundary method, we introduce a domain parameter tp to 
incorporate boundary conditions into the original governing equation. As men- 
tioned previously, the level sets of this domain parameter ip describe the diffuse 
boundaries and should satisfy n = V'?/'/|V'?/'|. On dQ, we impose a contact an- 
gle, 9, such that n- (V0/|V0|) = —cos 6', where V0/|V0| is the unit normal 
vector of the phase boundary (pointing to regions where = 1). We can thus 
derive the following equation for the boundary condition: 



Vt/-- V0 = -|V^|cose— . (20) 

This contact-angle boundary condition is similar to that suggested by Warren 
et al. [19] for contacting a sharp interface, in which a Dirac delta function 
replaces \Vip\. 

The chemical potential that drives the morphological evolution is defined by 
the variational derivative of the total free energy of the system: /i = 6F/6(f) = 
df/d(j) — e^V^0. We can apply the smoothed boundary formulation to the 
chemical potential by multiplying it by the domain parameter ip and applying 
the product rule to obtain: 

df e^ e^ df e^ elV^I r- 

where Eq. (20) was used in the third term. 

For a nonconserved order parameter in the phase field models, the evolution is 
governed by the Allen-Cahn equation [36] (also known as the time-dependent 
Ginzburg-Landau equation [37]), in which the order parameter evolves accord- 
ing to the local chemical potential variation: 

|._M,^-M(|-^V.(.V,)-fyV^cosA (22) 
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where M is the mobihty coefficient. 



For a conserved order parameter, the evolution of the order parameter is gov- 
erned by the Cahn-Hilhard equation, where the rate of the order parameter 
change is equal to the divergence of the its flux, which is proportional to the 
gradient of the chemical potential [Mf35] . d(j)/dt = V ■ (MV/x). The smoothed 
boundary formulation (derived in a similar manner as in the derivation in 



Section Q is given by i){d(p/dt) = V ■ {ipMVjj.) - V^ ■ (MV/x). Note that 
—MVfi = j is the flux of the conserved order parameter; thus, the second term 
represents the fluxes normal to the domain boundary (equivalent to Eq. (|5|). 
The smoothed boundary formulation of the Cahn-Hilliard equation is thus 
written as: 



d(j) _ 1 



V- 



^MV| 



^P 



V ■ (V^V0) 






2 f cos 9 



+ !^ J„. (23) 



where Jn = n ■ j. For a closed system, J„ is zero. Note that a small nonzero 
value should to be added to the domain parameter in the denominators in 



Eqs. (22) and (23) to avoid division by zero. 



4 Validation of the Presented Approach 



We herein demonstrate the validity and accuracy of the approach introduced 



in Sections 3.2 and 3.3 and the phase transformation with the presence of 
additional surface in Section 13.51 



4-1 ID Diffusion Equation 



First, we performed a ID simulation to demonstrate that the Neumann and 
Dirichlet boundary conditions were satisfied on two different sides of the do- 
main. Fick's second diffusion equation with given source and sink terms was 
solved within the domain defined by ■?/' = 1. The diffusion coefficient was set 
at 1, and the source and sink strengths were 0.02 and 0.01, respectively. On 
the right boundary of the diffusion domain, the gradient of C was set at -0.1, 
and on the left boundary, the value of C was set at 0.4. We selected the ID 
computational box for < x < 30 and used a hyperbolic tangent function for 
the continuous domain parameter ip: 



^ 



, /x-5\ , /x-25 
tanh — : — — tanh : — 



(24) 
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where ( is the coefficient for adjusting the interfacial thickness. The interfacial 
thickness is given approximately by ^o = 4.185^ where the interfacial region 
is defined by the range, 0.015 < ip < 0.985. The left and right interfaces are 
located at x = 5 and x = 25, respectively. We applied the smoothed boundary 
formulation, as in the derivation of Eq. (|8|, to reformulate the original diffusion 
equation, dC/dt = d^C/dx^ - C/0.01 + 0.02, to: 



dC 
~dt 



il) dx\ dx 



0.1 



dtjj 



1 



dx dx 



dx 



0.4 



Hv{ln 



dil) 



dx 



[1 - Hv{lm)] 



C 

(Toi 



(25) 



+ 0.02, 



where Hv{lm) is the Heaviside function used to specify the choice of the bound- 
ary condition and /^ = 15 is the midpoint of the diffusion domain. Therefore, 
the second and third terms only apply to the right and left interfaces, respec- 
tively. The initial concentration was C = everywhere in the computational 
box. A standard central finite difference scheme in space and an Euler explicit 
time scheme were employed in the simulations. 



Figure ^ shows the concentration proffies recorded at four different times (in 
solid blue lines). The domain parameter is plotted in the red line (the red 
circular markers indicate the position of grid points). The computational 
box was discretized to 1,200 grid points (Ax = 2.5 x 10~^), and ( was 
taken to be 2.86 x 10~^, such that the interfacial thickness is approximately 
^0 = 0.1197 = 4.79Ax. The parameters are given as Case lb in Table [l} On 
the right interface, it can be clearly observed that dC /dx = —0.1 at all times 
(except for a rapid change from dC/dx = to dC/dx = —0.1 in the very 
early transient period). In the early period, the concentration even took neg- 
ative values to satisfy the gradient boundary condition imposed at the right 
interface. In contrast, the concentration remained at 0.4 at the left interface 
during the entire diffusion process (except in the very early transient period, 
during which C changed from to 0.4). The analytical solution for the origi- 
nal sharp interface equation is also plotted for comparison, showing excellent 
agreement between the two methods. This result clearly demonstrates that 
both Neumann and Dirichlet boundary conditions are satisfied on the diffuse 
interfaces, and the smoothed boundary formulated equation reproduces the 
same result to the corresponding sharp interface version. 

To further analyze the effects of interfacial thickness and discretization resolu- 
tion on the smoothed boundary method, various simulations were conducted. 
In the first case, various interfacial thicknesses were selected, as in Case 1 
in Table ^ while the grid size was kept at Ax = 2.5 x 10"^. Figure p[a) 
shows the concentration distributions at t = 1,000 (nearly equilibrium) for 
C = 2.86 X 10-2 (Case lb) and C = 4.58 x 10"^ (Case If), for which the in- 
terfaces approximately span 4.79Ax and 76.7Ax, respectively. It is clear that 
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the calculated concentration deviates farther from the analytical solution when 
the interfacial thickness is greater (as shown in the derivation in Appendix [A)) . 
Figure |3tb) illustrates the relative errors during concentration evolution for 
various values of (. Here, the relative error is defined by the root-mean-square 
deviation (between the smoothed boundary result and the analytical solution) 
divided by the average analytical concentration. The results clearly show that 
the error increases as the interfacial thickness increases. As the concentration 
evolution approaches equilibrium, the errors also converge to their equilibrium 
values, as listed in Table [ij A scaling of the error to the interfacial thickness 
is observed as ( is varied from 2.86 x 10^^ to 4.58 x 10^^. In addition, the 
deviation between the smoothed boundary results and the analytical solution 
is much larger near the left boundary than near the right boundary, indicating 
the error associated with a Dirichlet boundary condition is larger than that 
with a Neumann boundary condition; see Figs. [3]^c)-(d). 

In the second case, we examined the effect of varying Aa; without changing the 
number of grid points across the interface. This was accomplished by selecting 
various grid sizes while maintaining the ratio of interfacial thickness to the 
grid size at 4.79. Results similar to Case 1 were obtained (see Case 2 in Table 
[I] and Figs. |4Fa)-(b)). Error increases with Ax. Since the resolution of the 
interface is unchanged (i.e., the number of points across the interface is fixed), 
it implies that the increased interfacial thickness is the dominant source of 
error. However, the errors are in general smaller than those in Case 1, which 
may be due to the fact that the interfacial thickness is effectively reduced when 
the resolution is decreased (the parts of the interfacial regions where ip is near 
and 1 are not resolved by large grid spacing). The same reason may explain 
the steep drop in the error for ( = 1.43 x 10~^ in Case la in Table [T| for which 
the rapid transition of ip is not properly resolved by the discretization. 

In the third case, we selected various grid sizes to examine the effect of the res- 
olution across the interface while maintaining the interfacial thickness (specif- 
ically, fixing ( at 5.73 x 10~^). The results show that error decreases when 
a larger grid size is selected; see Case 3 in Table [I] and Figs. |4|^c)-(d). This 
can be understood as follows. As we observed earlier, the smoothed bound- 
ary formulation reduces to the bulk partial differential equation far from the 
interface, where the gradient of the domain parameter vanishes. In the inter- 
facial region, the bulk term and the boundary term together set the boundary 
condition, as shown in Appendix |Xj In between, there is a region where the 
bulk equation is affected by the boundary term, which is small because the 
gradient is small, but not negligible. When the resolution is sufficiently low, 
the domain parameter in these regions take the bulk values and vanishing 
the boundary term and thus increasing the accuracy. As can be observed in 
Fig. |4Fd), the error behavior in such case is very different from other cases. In 
the specific example presented here, the discretized interface at the low res- 
olution of Ax = 0.2 is nearly a Heaviside step function, which yields smaller 
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error than the high-resolution cases. (When the resolution is high enough, the 
error is not affected by the resolution.) Therefore, in the ID case, an interface 
does not need to be fully resolved, and in fact the accuracy can be increased 
by not doing so. However, we found that numerical instability ensues when the 
resolution is further reduced. It has been determined that at least one point 
with an intermediate value between the two bulk values is required in order 
to achieve numerical stability. 

Above argument applies to only ID case or when interfaces are very flat in 
multiple dimensions. In one dimension, the curvature of the interface is zero 
(i.e., the interface is flat), and therefore the effect of curvature can be neglected 
and a good resolution across the interface (which provides the smoothness of 
the curved interface) is not required. This is not the case when sufficiently 
large curvature is present, and thus smoothed interface with about three grid 
points are required to obtain accurate results for 2D or 3D calculations. 

In the fourth case, we varied the value of v (the small value added to the 
denominators to avoid division by zero) from 1 x 10~^ to 1 x 10~^^, while the 
grid size and interfacial thickness were maintained at Ax = 2.5 x 10~^ and 
^0 = 4.79Aa;, respectively. In practice, a smaller v would lead to a less stable 
numerical implementation because the values of l/ip or l/ip"^ become much 
larger, which requires a much smaller time step size. The results show that 
the error quickly converges to a small value when v is smaller than 1 x 10"^; 
see Case 4 in Table [l] This suggests that once v is small enough to yield 
converged results, further reduction is unnecessary and should be avoided so 
that a larger time step can be employed. 

In summarizing the above ID test simulations, we found that the interfacial 
thickness is the dominant source of error. The errors are less sensitive to the 
resolution of the finite-differencing discretization (selection of Ax) and the 
parameter for singularity control (selection of v). When the diffuse interface is 
properly resolved, the error scales with the interfacial thickness. Moreover, in 
general, the error that results when a Dirichlet boundary condition is imposed 
is larger and more sensitive to the interfacial thickness than when a Neumann 
boundary condition is imposed. This behavior can be understood from the 
results of analysis in Appendix lAl where the scaling of the errors can be found 



in Eqs. (A. 3) and (A.5). 



4^.2 Surface Diffusion and Bulk Diffusion in a Cylinder 



To further demonstrate the validity of the smoothed boundary method, we 
applied the method to simulate oxygen-vacancy diffusion in a cylinder, for 
which a cylindrical coordinate grid system was used. We solved the coupled 
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surface-bulk diffusion problem using both the smoothed boundary and the 
original sharp interface formulations in the same grid system for comparison. 
For the smoothed boundary method, we used a hyperbolic tangent function, 
ilj{r,z) = {1 — tanh [(i? — r)]/C}/2, of the continuous domain parameter to 
define a cylinder, where r is the radial position, z is the axial position, and 
R is the cylinder radius. Therefore, the cylinder surface {ip = 0.5) where 
surface reaction and surface diffusion occur, is located at r = R, the solid 
region ('0 = 1) for bulk diffusion is defined at r < i?, and the environment 
{■ip = 0) is defined at r > _R. We selected the cylinder radius R and the 
cylinder axial length to be 1 and 12, respectively. The grid sizes were selected 
to be Ar = 1.76 x 10~^ and Az = 4 x 10~^, such that the cylinder contains 
57 and 300 grid points in the radial and axial directions, respectively. (The 
computational box is larger than the cylinder in the radial direction, and 
contains 75 and 300 grid points in the radial and axial directions, respectively.) 
The interfacial thickness was selected to be 4.26Ar by setting ( = 1.0182Ar. 



Equation ( 12 ) was solved using a standard central finite difference scheme in 
cylindrical coordinate system and an Euler explicit time scheme; see Appendix 
[B] for the discretization scheme. The parameters for the diffusion equation 
were selected to be Dy, = 1, Dg = 10, and L = 0. The boundary conditions 
at the computation box boundary were set at C = 1 at z=0 and C = at 
z = 12, with no gradient on the remaining two sides. For comparison, the 



original sharp interface equations, Eqs. (^ and (10), were solved using the 
same discretization scheme with the same grid system and resolution. For this 
case, the surface-reaction-diffusion boundary condition, Eq. ([9]), was explicitly 
imposed at the 57th grid points in the radial direction. The concentration 
evolution is implemented as follows. First, the surface concentration is updated 
by Eq. ^ according to the normal flux at the cylinder surface calculated 



from the normal gradient of the surface concentration obtained from Eq. ( 10 ) 
Next, the normal surface flux is calculated using Eq. ^ with the updated 
surface concentration. The cylinder concentration is then evolved according to 



Eq. (10) with the normal flux boundary condition. This procedure is repeated 
within the Euler explicit time scheme for the concentration evolution. For 
the smoothed boundary formulation, we simply solve a single equation that 



automatically includes coupled bulk and surface diffusion, Eq. (12). 



Figures |5](a) and (b) show the steady-state concentration profiles of the sharp 
interface version for k = 2.1 and k = 50, respectively. The concentration 
decays along the axial direction according to boundary values prescribed at the 
box boundaries. The diffusion front bends because of the surface reaction, such 
that the concentration is lower near the cylinder surface. Shown in Figs. [51(c) 
and (d) are the corresponding smoothed boundary results. For clarity, only 
the concentration in the region of < 2; < 6R is presented. The results from 
the two methods are in excellent agreement, clearly demonstrating the utility 
and validity of the smoothed boundary method for incorporating two sharp- 
interface equations into one smoothed boundary equation. 
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To further examine the effect of interfacial thickness, we included two other ra- 
dial grid sizes in the simulations, i.e., Ar = 3.49 x 10~^ and Ar = 6.86 x 10~^, 
such that the cylinders contain 29 and 15 radial grid points, respectively. By 
selecting the radial grid sizes in this way, each radial grid point in a lower- 
level resolution (thicker interface) overlaps with every other grid point in a 
higher- level resolution (thiner interface). The diffuse interface is maintained 
to span 4.26Ar, and the axial grid size is kept at Az = 4 x 10~^. Hereafter, we 
refer the three interfacial thicknesses to as the "thin" (^o = 0.075), "medium" 
(^0 = 0.149) and "thick" (^o = 0.292) interface cases. Shown in Figs, ge) and 
(f) are the steady-state concentration profiles for the thick-interface results, 
corresponding to the cases in Figs. [5ta) and (b). The results are still in rea- 
sonably good agreement with the original sharp interface results, even though 
the interfacial thickness is approximately 29.2% of the cylinder radius. 

The relative errors of the thin, medium, and thick interface smoothed bound- 
ary results are plotted in Fig. |6j The relative errors are calculated by dividing 
the differences between the smoothed boundary and sharp interface results 
by the average concentration of the sharp interface results. The average con- 
centration is calculated for the active region between the plane at 2r = and 
the plane on which the maximum concentration is 0.01. Note that only the 
errors within in the cylinder defined by V' > 0.5 are considered. For the thick- 
interface smoothed boundary result, the maximum local errors for k = 2.1 and 
K = 50 are approximately 5 x 10~^ and 0.075 (see Figs. |6](e) and (f)), whereas 
the average relative errors are 1.81 x 10"'^ and 2.90 x 10~^, respectively; see 
Table |2] The average relative errors, denoted by e in Table [2| are calculated by 
dividing the root-mean-square deviation between the smoothed boundary and 
the sharp interface results by the average sharp interface concentration in the 
cylinder. The root-mean-square deviation and average concentration are cal- 
culated in the cylindrical coordinate system. As expected, the error increases 
as the interface becomes thicker (i.e., as ^o increases). However, in contrast 
to the ID simple diffusion test in Section |4.1[ the behavior of the error is in- 
consistent across the parameter sets; see Table |2j This relatively complicated 
error behavior may originate from the coupling of the bulk and surface diffu- 
sion equations. In addition to the effect of the interfacial thickness, the error 
also increases with a larger reaction coefficient k, which may be explained by 



the increase of the scaling coefficient for the error (/iq defined above Eq. (A. 2) 
in Appendix lAl) when the given boundary value is larger. In addition, the 
gradient of the concentration near the boundary increases in magnitude with 
increasing boundary condition value, which can lead to a larger error. 

One interesting phenomenon is observed in the high k results (k = 50 and 
100). Although the errors in the bulk greatly increase with the interfacial 
thickness, the errors at the surface remain small; see Fig.|6](f) and Table[2] This 
indicates that the error originates from the boundary condition affecting the 
bulk solution, rather than from an increased error in the boundary condition 
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value. The thicker interface thus leads to a larger bulk region that is affected by 
the boundary condition. Therefore, we compare the error associated with the 
bulk region and with the boundary condition. Here, the bulk errors, denoted 
by eb in Table [2} are calculated by the same method as the average relative 
error but exclude the grid points on the nominal cylinder surface {r = R). The 
surface errors, denoted by Cg in Table |2| are calculated by the same method but 
with only the grid points at the nominal cylinder surface. In the case where 
K = 100, the error at the surface even decreases with interfacial thickness. 



4-3 Contact-Angle Boundary Condition 



We performed simple 2D simulations to validate the smoothed boundary for- 
mulation for the contact-angle boundary condition at the three-phase bound- 



ary. Equations (22) and (23) were tested for nonconserved and conserved order 
parameters, respectively. The equations were solved using the central finite dif- 
ference scheme and the Euler explicit time scheme. The computational box 
sizes are L^ = 100 and Ly = 100, and the parameters used are Ax = 1 and 
M = 1. A simple common double- well function was selected for the bulk free 
energy functional, /(0) = w(j)'^{l — 0)^, such that the steady-state phase field 
order parameter profile is determined by = {1 — tanh [(A/wx)/(A/2e)]}/2, 
where x is the coordinate variable indicating the distance to the phase bound- 
ary, and the characteristic thickness of the diffuse interface is determined by 
6^ = eJ2/w. By setting e = Jl/w, the characteristic thickness is controlled by 

S^ = a/2/w, and the phase field interfacial energy is maintained at a constant 
value, 70 = e\/2w/6 = v^/6. A horizontal diffuse-interface flat substrate sur- 
face is defined by the hyperbolic tangent function ip = {1 -|-tanh [{y — 30) /C]}, 
such that ip = 0.5 is at y = 30, and ip gradually transitions from to 1 from 
below to above the substrate surface. Here, we have two diffuse interfaces: one 
for the phase field order parameter and the other for the smoothed boundary 
domain parameter. Both thicknesses can affect the accuracy when imposing 
contact-angle boundary conditions. To verify the contact angle boundary con- 
ditions, various combinations of substrate surface thickness and phase bound- 
ary thickness were selected by adjusting the values of ( and w. The initial 
phase boundary was placed vertically in the middle of the domain (x = 50), 
with Phase 1 (0 = 1) and Phase (0 = 0) on the left and right halves, re- 
spectively. On the computational box boundaries, the normal gradients of the 
phase field order parameter were set at zero: d(j)/dx = at x = and 100, 
and dcp/dy = at ?/ = and 100, which can be interpreted as the no-fiux 
boundary conditions. 



In the first set of simulations, we evolved Eq. (22) for a nonconserved order 
parameter with a 60° contact angle. The result clearly shows a 60° contact 
angle at the three-phase boundary, as specified; see Fig.JTJ^a). The angle can be 



measured at the intersection between the two contours oiip = 0.5 and = 0.5, 
as shown in Figs. Wih) and (c). The 60° angle is maintained during the entire 
evolution, except for the very early transient period, when the contact angle 
changes from the initial 90° angle to the prescribed 60° angle. Because of 
the contact-angle boundary condition, the initially flat phase boundary bends 
and creates a negative curvature in Phase 1. As a result, the phase boundary 
moves toward Phase 0. Once the phase boundary evolves to a circular arc 
with a uniform curvature everywhere (other than regions in contact with the 
substrate), it moves at a uniform constant speed in a steady-state motion, and 
eventually only Phase 1 remains in the system. 



In the second set of simulations, we evolved Eq. (23) for a conserved order pa- 
rameter in a closed system with J„ = and a 120° contact angle. As expected, 
the phase boundary intersects the substrate surface at a 120° contact angle; 
see Fig. [7Fd)-(f). In contrast to the Allen-Cahn-type dynamics, because of the 
conservation of the order parameter, the phase boundary near the substrate 
moves toward the left, whereas the phase boundary away from the substrate 
moves in the opposite direction. As a result, the phase boundary deforms into 
a curved shape. When the system reaches its equilibrium state, the phase 
boundary forms a circular arc with a uniform curvature everywhere (except 
where the phase is in contact with the substrate), such that the total surface 
energy is minimized; see Fig. uld) for t = 3.0 x 10^. 



Table [3] lists the average values of cos 6' calculated by {Vip ■ V0)/(|V'?/'||V0|) 
at the grid points within the three-phase boundary region defined by 0.1 < 
ip < 0.9 and 0.1 < cf) < 0.9 in the steady state (i.e., when the phase bound- 
ary becomes a circular arc). These results again clearly show that the error 
(for a given phase boundary thickness) increases as the interfacial thickness 
increases. This can be understood based on the analysis of the ID test results 
in Section Wa\ since the substrate surface is assumed to be flat in this test. If 
the substrate interface is curved, the resolution of the interface will have more 
influence on the error. 

In contrast to the effect of domain boundary thicknesses, the error is relatively 
insensitive to the phase boundary thickness once the phase boundaries are 
properly resolved; see cases with 5^ >= 1.4142 in Table [3] However, when 
the phase boundary is too thin, the error tends to increase because of the 
loss of resolution in the phase- field-order-parameter gradient; see cases with 
5^ = 1.0607 in Table [31 In general, the results demonstrate that the contact- 
angle boundary condition is well imposed using the presented method. Even 
when the domain boundary thickness is as large as 16.74 grid spacings with 
C = 4, the contact angle only deviates less than 2° from the imposed values 
as long as the phase boundary is properly resolved. Here, the error is sensitive 
to the resolution of the interface because the phase boundary is curved unlike 
the substrate surface. 
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In addition to the contact angle, a no-flux boundary condition for a conserved 
order parameter is implicitly imposed at the substrate surface. The error asso- 
ciated with such a boundary condition was evaluated by examining the overall 
change in the value of the total order parameter. The conservation of the order 
parameter was met within a numerical error (well below 1% in most cases) in 
these validation simulations; see Table [3] 



5 Applications 



Although the details of the scientific calculations performed applying these 
methods to problems in materials science will be published elsewhere, it is 
worth presenting some of the results herein to demonstrate the potential of 
the method. 



5.1 Oxygen- Vacancy Diffusion in SOFC Cathode 



The first example is ionic transport through a complex microstructure. Here, 
ion diffusion is driven by a sinusoidal voltage perturbation. For the steady- 
state solution, the time dependence of the form exp(ici;t), where u is the an- 
gular frequency and i = a/— 1, can be removed as in the equation derived by 
Lu et al. [lOj. For this case, the smoothed boundary formulated equation is 



obtained from Eq. (12) to: 



V ■ (tpD^VC) - \V^\{kC - D.VlC) = i^uC, (26) 

where C is the complex concentration amplitude, which consists of real and 
imaginary parts and includes the amplitude of the concentration wave and 
the phase shift. Note that the surface accumulation term is ignored {L = 0) 
here because its magnitude is usually very small in comparison with the bulk 
concentration [40j. This equation can be solved by an alternating direction 
line relaxation (ADLR) method in a second-order central-difference scheme in 
space; see Appendix [B] for the numerical implementation. 

In this work, we adopted an experimentally reconstructed complex microstruc- 
ture, the porous ceramic cathode and nonporous ceramic electrolyte of an 
SOFC, as the input geometry. The microstructure data is stored as a 3D array 
consisting of 321 x 261 x 297 voxels that indicate the electrolyte (gadolinia- 
doped ceria: GDC), cathode (lanthanum strontium chromite: LSC) and pore 
phases by different values. To emphasize the convenience of image-based smoothed 
boundary simulations, we treat the center of each voxel as the location of the 
grid points in the calculation without further enhancement of the resolution 
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from our initial reconstructed microstructure. For very high-accuracy scientific 
calculations, one can easily enhance the resolution by refining the grid sizes. 
To smooth the voxelated, discontinuous data, we first employed a level set 
distance function method |t42j to determine the distances between grid points 
and the solid-pore interface, and then computed the hyperbolic tangent of the 
distance function to obtain the domain parameter profile; see Appendix [E] for 
details. 

For simulations of the concentration distribution in the porous cathode, the 
regions containing nonporous electrolyte are excluded, such that the compu- 
tational box only consists of 321 x 176 x 297 grid points. The grid spacing 
is set at Ax = 6.285 x 10^^. The boundary conditions along the main dif- 
fusion direction (the y-axis) on the computational box are Re(C) = 1 and 
Im(C') = at y = 0, and Re(C') = and lm{dC/dy) = at y = 11.062. 
The boundary conditions on the remaining four sides are zero-gradient for 
both the real and imaginary parts. As a demonstration of the method, the 
length scale and physical material properties are nondimensionalized. Figure 
[8] shows the steady-state concentrations for the cases in which surface diffusion 
is excluded (k = 0.1 and Dg = 0) and included (k = 2.1 and D^ = 10) with 
-Db = 1 and direct current (DC) loading {u = 0). In these cases, the imagi- 
nary part vanishes, and the solution of the real part is equivalent to that of 



a homogeneous Helmholtz-like equation with the right-hand side of Eq. (26) 
equal to zero. As shown in Fig. |8} the concentration decays from 1 to along 
the j/-axis over the complex cathode microstructure to satisfy the boundary 
conditions imposed on the box boundaries and at the cathode-pore interfaces. 
The utilization lengths (i.e., the length over which the cathode material is 
active) of the two cases are similar, as predicted by Lu et al. |10] for a cathode 
with simplified cylindrical geometry, in which the effective diffusivities under 
DC loading with and without surface diffusion are found to be similar for the 
parameters given above. However, a slight difference in the concentration dis- 
tributions of the two cases can be observed. Because of the faster transport 
path along the surface, the diffusion front with surface diffusion (Fig. [sFb)) is 
more planar compared with that without surface diffusion (Fig. |8](a)). 

Figure [9] shows the real and imaginary parts of the steady-state concentration 
amplitude for the cases in which Db = 1, k = 2.1 and Dg = 10, with alternating 
current (AC) loading of the angular frequencies oi u = 1.5 and 51.5. The 
boundary conditions on the computational box are the same as in the DC 
loading case above. In the low frequency case (Figs. [OJ^a) and (d)), the real 
part of the concentration, which represents the amplitude of the concentration 
wave decays and forms a planar diffusion front within the utilization length, 
where the material is active {0 < y < 5). Additionally, a negative value of 
the imaginary part occurs in the regions where the real part decays because 
of the phase shift resulting from the delayed response. The magnitude of the 
imaginary part then decays back to zero toward the inactive region. In the 
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high frequency case, the enhancement of concentration along the surface is 
observed due to surface diffusion; this is evident in Figs. loFb) and (c), which 
show larger values of the real part of the concentration amplitude within the 
utilization length {0 < y < 2.5). The real part of the concentration amplitude 
quickly decays from the surface into the bulk. In contrast to the enhanced 
real part at the cathode surface, the magnitude of the imaginary part is small 
near the surface because surface diffusion reduces the response time and the 
phase change is thus decreased. In an analogy to the low-frequency response, 
a negative imaginary part occurs in the region where the real part decays. 
The magnitude of the imaginary part decays toward the inactive region. This 
behavior can be more clearly discerned in the magnified views in Figs. |9](c) 
and (f). 

The smoothed boundary method can also be used to impose Dirichlet bound- 
ary conditions on irregular surfaces. For example, if the ionic diffusivity in 
the electrolyte is assumed to be much larger than that in the cathode, the 
concentration in the electrolyte will be nearly uniform. To simulate this sce- 
nario, we impose a fixed concentration at the electrolyte-cathode contacting 
surface as the boundary condition. We used the experimentally reconstructed 
321 X 261 X 297 array that contains a porous cathode and a nonporous elec- 
trolyte as our input geometry. The voxelated data were smoothed to the hy- 
perbolic tangent domain parameter profile by the level set distance function 
method mentioned in Appendix |E} Here, three domain parameters are em- 
ployed to define the three regions: electrolyte (ipi), cathode (V'2), and pore 
{ip^ = 1 — ipi — ip2)- The smoothed boundary formulated governing equation 



is obtained by modifying Eq. (12) to: 



dC V ■ (^2/^bVC) |V^2 



dt 4)2 1p2 

^b 



kC - DS7tC 



Wn- 



^l 



v^2-v(^2C7)-|v^2rs^ 



(27) 



W, 



D, 



where the weighting factors are given by Wn = [|V'?/'2||VV'3|/(|VV'i||V^2| + 
|VV'2||VV'3| + |V7/;3||V^i|)]'^andiyz) = [|V^i||VV'2|/(|VV'i||VV'2| + |V^2||VV'3|+ 
iV^sllV^il)]^, such that the Neumann boundary condition (surface reaction 
and surface diffusion) is imposed only at the cathode-pore interface (| VV'2| I Vt/'sI 7^ 
0), and the Dirichlet boundary condition (a prescribed concentration value) is 
imposed only at the electrolyte-cathode interface (|V'?/'i||V'?/'2| 7^ 0). The expo- 
nent P determines the transition profiles from the Neumann to the Dirichlet 
boundary conditions in the regions of three-phase boundaries. We selected 
(3 = 0.8 for this numerical simulation. On the computational box boundaries, 
we set C = at y = 16.404 and the zero-gradient boundary condition for the 
remaining five sides. 

The same material parameters used in the cases of Fig. [8] were selected. Fig- 
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ures 10 (a) and (b) illustrate the reconstructed SOFC complex microstructure 
and irregular surfaces defined by the values and gradients of the domain pa- 
rameters, respectively. Figures [lotc) and (d) show the simulation results of the 
steady-state oxygen-vacancy concentration distributions with a fixed value of 
C = 1 imposed at the cathode (LSC)-electrolyte (GDC) interfaces. The con- 
centration distribution is very different from the ones shown in Fig. |8] because 
a larger portion of lateral diffusion occurs in the x and z directions, which re- 
sults from the smaller contacting areas (compared to the cross-sectional area 
of LSC on the x-z plane in Fig. pi where diffusion is mainly in the y direction). 
As a result, the concentration drops rapidly within a short distance from the 
contacting areas, making the utilization length of the cathode material shorter 
and uneven. 



5.2 Kirkendall- Effect Diffusion with a Moving Boundary 



5.2. 1 Kirkendall-Effect-Induced Deformation Modeled by Navier- Stokes- Cahn- 
Hilliard Equations 



The third application demonstrates the smoothed boundary method's broad 
applicability by applying it to the coupled Navier-Stokes-Cahn-Hilliard equa- 
tions |43I44|45|46|I47|I48] . This particular formulation aims to solve diffusion 
problems with the Kirkendall effect with efficient and abundant vacancy sources 
and sinks in the bulk of a solid |1^50|51II52|I53] . In this case, the solid ex- 
periences deformation because of vacancy generation and elimination. The 
Navier-Stokes-Cahn-Hilliard equations are coupled with the smoothed bound- 



ary formulation of the diffusion equation derived in Section 3^ as a model 
of plastic deformation because of volume expansion and contraction resulting 
from vacancy fiow. 



When the diffusing species of a binary substitutional alloy have different mo- 
bilities, the diffusion fluxes of the two species are unbalanced, creating a net 
vacancy flux toward the side containing the fast diffuser. Here, we denote the 
quantities associated with the slow diffuser, fast diffuser and vacancy by the 
subscripts A, B, and V, respectively. Because of the accommodation/supply 
of excess/depleted vacancies, the sohd locally expands/shrinks [M|55p56f57f58j 
when maintaining the vacancy mole fraction at its thermal-equilibrium value. 
We treat the solid as a very viscous fluid j5^6Uf61IIB^63j with a much larger 
viscosity than that of the surrounding environment. In this case, we solve the 
Navier-Stokes-Cahn-Hilliard equations to update the shape of the material as 
follows 



-VP + V -T] 



Vv + (Vv) 



^[i'^l^i^^^ 



(28a) 



23 



V ■ V = (?v., (28b) 

|:-v.V^ = MV^(|^-e^V^V^), (28c) 

where P is the effective pressure, r] is the viscosity, v is the velocity vector, d is 
the number of dimensions, the superscript T denotes the transpose, C^ is the 
Cahn number reflecting the capillary force compared to the pressure gradient, 
gv is the vacancy generation rate per unit volume, and if) is the domain param- 
eter indicating the solid phase for diffusion. One great advantage in employing 
a phase-field type equation is that it automatically maintains the profile of the 
domain parameter, ip, in the form of a hyperbolic tangent function because 
it is the equilibrium solution for the phase field equation (Eq. ( |28c )). Note 



that here we ignore the inertial force in the Navier-Stokes equation to ob- 



tain Eq. (28a) because the deformation is assumed to be a quasi-steady-state 
process. The vacancy generation rate that results in the local volume change 
(dilatational strain) is given by gv = —[V ■ {Dvb^Xb)]/[pi{^ — Xy)], where 
Xb is the mole fraction of the fast diffuser, Xy is the thermal-equilibrium 
vacancy mole fraction (which is assumed to be maintained throughout the 
solid in this model), Dys is the diffusivity for vacancy flux associated with 
VXb, and pi is the lattice site density of the solid. Here, Xy is taken to be 
1.6 X 10~^. The evolution of the fast diffuser mole fraction is governed by the 
advective Pick's diffusion equation, written as: 



dt 



V . VXb = V ■ {DIj^VXb) - Xsgy. (29) 



where -D]^^ is the diffusivity for the fast diffuser flux associated with VX^, 
and the advective term accounts for the lattice shift because of volume change. 
Because diffusing atoms cannot depart from the solid region, a no-flux bound- 
ary condition is imposed at the solid surface. Thus, the smoothed boundary 



formulation of Eq. (29) is written as: 



^-v.V.Y.^ ^-W^««^^--) -A-.,.. (30) 

at ip 



As the concentration evolves, the shape of the solid is also updated by Eq. ( 28c ) 



and by iteratively solving Eqs. (28a) and (|28b) through the application of a 



projection method [65|66] : see Appendix [F] for the numerical implementation 

The slow and fast diffusers are initially placed in the left and right halves of 
the solid, respectively. We use their theoretically calculated diffusivities for 
this simulation [E7|68|69ll7n] . Here, we calculate the slow diffuser atomic hop 
frequency based on the material parameters of aluminum at 600 K, and set 
the fast diffuser atomic hop frequency four times larger than that of the slow 
diffuser. Figure [IT] shows snapshots of the mole fraction profiles (left column) 
and velocity fields (right column) from a 2D simulation. As the fast diffuser 
diffuses from the right to the left side, the vacancy elimination and generation 
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cause contraction and expansion on the right and left sides, respectively. As a 
result, the initially rectangular slab deforms into a bottle-shaped object. 



5.2.2 Kirkendall Void Growth with Localized Vacancy Sources 

In another scenario in which the vacancy diffusion length is comparable to 
or smaller than the distance between vacancy sources and sinks, the explicit 
vacancy diffusion process must be considered [^70II71|I72] . In this case, va- 
cancies diffuse in the same manner as the atomic species. In the bulk of a solid 
devoid of vacancy sources/sinks, the concentration evolutions are governed by: 



dX 



V 



V ■ {Dvv^Xv + Dvb^Xb), (31a) 

V ■ {DbvVXv + DI^VXb). (31b) 



dt 

dt 

Because the solid surfaces are very efficient vacancy sources/sinks [6^f72] . we 
impose the thermal-equilibrium vacancy mole fraction at the solid surfaces 



as the Dirichlet boundary condition for solving Eq. (31). In this case, the 



smoothed boundary formulation of Eq. (31) is given by: 



^ = Iv . [^{DvvVXv + DvbVXb)] - ^, (32a) 



^ = ^V ■ [^{DBvVXy + DI^VXb)] + Y^^2^ (32b) 

where K = DyylVip ■ Vi^ipXy) — \Vip\'^Xy]. Because the vacancy generation 
and elimination in this scenario only occurs on the solid surfaces, internal 
volume change in the bulk is not considered. Therefore, instead of using a 
plastic deformation model as in the previous case, we adopt a typical Cahn- 
Hilliard type dynamics to model the shape change: 






MV^^-eV^^ +^.-^, (33) 



where Jy = —{DyyVXy + DvbVXb) is the vacancy flux, and the last term 
represents the normal velocity of the solid surfaces because of vacancy injection 
into or ejection from the solid. 

An example of the results obtained using this approach is the growth of a void 
in a rod [72|73|74|I75] . The above equations were solved using a central differ- 



ence scheme in space and an implicit time scheme (see Appendix G). The fast 
diffuser was initially placed in the central region while the slow diffuser filled 
the outer region. A void was initially placed off-centered in the fast-diffuser re- 
gion, where a ID study found to be the likely nucleation site [72]. The vacancy 



mole fractions were fixed at the void and cylinder surfaces. Figure 12 shows 
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snapshots of the fast diffuser mole fraction profile (normalized to the lattice 
density) and the vacancy mole fraction profile (normalized to its equilibrium 
value). As the fast diffuser diffuses outward, vacancies diffuse inward from the 
rod surface to the void surface, causing vacancy concentration enhancement 
and depletion in the central and outer regions, respectively. To maintain the 
equilibrium vacancy mole fraction at the rod and void surfaces, vacancies are 
injected and ejected at those surfaces. As a result, the rod radius increases, 
and the void grows. Similar dynamics were examined using a sharp interface 
approach [72], but this new method provides the flexibility in geometry to 
examine cases where a void initially forms off-centered. 



5.3 Thermal Stress 



Solid oxide fuel cells (SOFCs) usually operate at temperatures near 500° — 
1, 000°C. Evaluating the thermal stress resulting from the differences in ther- 
mal expansion and elastic moduli is important for analyzing mechanical fail- 



ure. We expand the generalized mechanical equilibrium equation, Eq. (17), for 
a linear, elastic and isotropic solid. (Note that the derivation for the mechan- 
ical equilibrium equation is general and is not limited to isotropic solids. We 
selected an isotropic model because of the lack of available crystallographic 
information among the experimental data.) The equation is discretized in a 
central finite difference scheme and numerically solved by an ADLR solver; 
see Appendix [H] for details. 

The thermal expansion rates of the ceramic electrolyte (GDC) and cathode 
(LSC) are taken to be 12.3 x 10"^ K"^ [76] and 10.6 x 10"^ K"^ [77], such 
that the thermal expansions at operation temperature are 0.0123 and 0.0106, 
respectively. (Here, we have assumed arbitrarily that the composite material 
is relaxed at a reference temperature, and assumed an operation temperature 
of 1000° above the reference temperature.) We chose the elastic constants of 
GDC to be isotropic (An — A12 = 2A44), and the values are Af/^*" = 375.94 
GPa, Af/^^ = 188.54 GPa and Af/^^ = 93.70 GPa, calculated from a Young's 
modulus of 250 GPa and a Poisson's ratio of 0.334 [751179] : see Appendix [h] 
The LSC phase is softer than the GDC phase, and its elastic constant is also 
assumed to be isotropic. The values are selected to be Aff*" = 269.23 GPa, 
\\l^ = 115.38 GPa and A^f^ = 76.29 GPa, based on a Young's modulus 



of 200 GPa and a Poisson's ratio of 0.3 [77]. As in Section 5.1, we again use 
domain parameters to indicate the GDC phase {ipi = 1 inside the GDC and 
t/'i = outside the GDC) and the LSC phase {4'2 = 1 inside the LSC and 
'?/'2 = outside the LSC). The entire solid phase is then represented by the 
sum of the two phases, ip = ipi + ip2 = '^^ The body force term and elastic 
constant tensor are replaced by an interpolated, spatially dependent thermal 
expansion and elastic constant tensor according to the domain parameters; 
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see Appendix [Hj The solid surface is assumed to be traction-free, Ni = 0. 



In this simulation, we selected the same computational box as in the case of 
Fig. [Io|(a), containing 321 x 261 x 297 grid points in the x, y, and z directions, 
respectively. Each grid point represents a voxel in the experimentally obtained 
microstructure. The yellow color indicates the LSC phase, and the semitrans- 



parent cyan color indicates the GDC phase; see Fig. 10 a). The grid spacing is 
Ax = 25 nm, such that the computational box spans 8.025 x 6.524 x 7.425 /im^. 
We assumed a rigid computational box with frictionless boundaries on the six 
sides, which means that u = dv/dx = dw/dx = on the two y-z planes, 
V = du/dy = dw/dy = on the two x-z planes, and w = du/dz = dv/dz = 
on the two x-y planes of the computational box boundaries, where u, v and 
w are the displacements along the x, y and z axes, respectively. While this 
set of boundary conditions is not realistic for SOFC material environment, we 
chose it for the demonstration purpose in order to avoid overlaps with a future 
publication of physically based SOFC simulations. 



Shown in Fig. |13[a) are the calculated mean stress distributions resulting 
from thermal expansion in a confined sample. The mean stress is defined 
by: (Tjn = {cxx + CFyy + o"2z)/3, whcrc the stress components are calculated 
according to the method provided in Appendix |Hj Here, we choose mean 
stress to illustrate the effective pressure in the solid. A negative mean stress 
indicates that the region is under compression. Despite a complicated stress 
distribution observed because of the complex geometry, the overall magnitude 
of the mean stress is roughly between 2 and 4 GPa, which can be roughly 
estimated by the product of Young's modulus and the thermal expansion 
with an enhancement resulting from the porosity of the solid. Additionally, an 
overall larger stress in the GDC phase is observed, reflecting the GDC phase 
is harder than the LSC phase. Figure [IsFb) shows the mean stress in the LSC 
phase and Fig. [Istc) shows the mean stress on the GDC surface after rotating 
the volume 180° around the 2;-axis. Three types of stress enhancements can 
be observed in the simulation result. At the cathode-electrolyte contacting 
surfaces, stress is enhanced because of the mismatch of thermal expansion and 
elastic constants between the two materials; see the red arrows in Figs. fl3tb) 
and (c). The second is the concentrated stress observed at the grooves on 
the electrolyte surface (not contacting the cathode), as shown by the white 



arrows in Fig. 13 c). The third type is the stress concentration effect at the 
bottlenecks in the cathode phase, where the stresses are roughly larger by a 
factor of three to four compared to the overall value, as shown by the green 
arrows in Fig. [l3Fb). The simulation results demonstrate that the smoothed 
boundary method can properly capture the linear elasticity and the geometric 
effects of the system based on a diffuse-interface defined geometry. 
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5-4 Phase Transformations in the Presence of a Foreign Surface 



The AUen-Cahn equation describes the dynamics of a nonconserved order pa- 
rameter, which can be taken as a model for the ordering of magnetic moments 
[38] and diffusionless phase transformations that involve only changes in crys- 
talline order [38] . This equation can also be used as a model for evaporation- 
condensation dynamics [38]l39] . Here, we use the AUen-Cahn equation to exam- 
ine the evaporation of a droplet on a rough surface. The domain parameter was 



given a ripple-like feature, as shown in Fig. 14, having a hyperbolic-tangent-like 
profile continuously transitioning through the substrate surface {i/j = 1 above 
the surface, and ip = below the surface). The droplet phase was placed on 
top of the boundary, and its shape was evolved by the smoothed boundary 



formulation of the Allen-Cahn equation, Eq. (22), using the standard central 
difference scheme in space and an Euler explicit scheme in time. The simula- 
tion was performed in two dimensions, using the parameters Ax = 1, M = 1 
and e = 1, with a domain size of L^ = 100 and Ly = 100. The contact angle 
was set at 135°, and a zero-gradient boundary condition of is set at the 
computational box boundaries. 



The evolution of the droplet surface as it evaporates is illustrated in Fig. 14 a) 
as a contour (0 = 0.5) plotted at equal intervals of 270 dimensionless time 
units. The color change from blue to red indicates various times from the ini- 
tial to the final stages, respectively. As the surface evolves, it is clear that 



the contact angle is maintained, as shown in Fig. 14 b). The dynamics of 
the motion of the three-phase boundary are interesting in that the velocity 
changes depending on the angle of the surface (with respect to the horizontal 
axis), which can be inferred from the change in the density of the contours. 
Because the interfacial energy is assumed to be constant, the droplet would 
prefer to have a circular cap shape. However, the contact angle imposes an- 
other constraint at the three-phase boundary. When the orientation of the 
surface is such that both of these conditions are nearly met, the motion of the 
three-phase boundary is slow as the droplet evaporates. When the orientation 
becomes such that the shape of the droplet near the three-phase boundary 
must be deformed (compared to the circular cap), the three-phase boundary 
moves very quickly, which leads to an unsteady motion of the three-phase 
boundary. In contrast, at the top of the droplet far from the substrate, the 
curvature is barely affected by the angle of the substrate surface; thus, the 
phase boundary there moves at a speed inversely proportional to the radius. 
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5.5 Motion of a Droplet Due to Unbalanced Surface Tensions 



In another example application, we modeled a self-propelled droplet. Here, 
two different contact-angle boundary conditions are imposed on the right and 
left sides of the droplet placed on a flat surface. The smoothed boundary for- 



mulation of the Cahn-Hilliard equation, Eq. (23), is used with J„ = in this 
simulation. The domain sizes are L^ = 240 and Ly = 60. The parameters and 



computational box boundary condition are the same as in Section |5.4[ The 
contact angle on the right side of the droplet is set to 45° and that on the left 
side to 60° by imposing position-dependent boundary conditions. Note that 
this setup is equivalent to the situation in which the substrate-environment, 
droplet-substrate and droplet-environment surface energies satisfy the condi- 
tions of Young's equation: 

7se — 7sd = 7de COS 60° for the left side, (34a) 

7se — 7sd = 7de COS 45° for the right side, (34b) 

where 7se, 7sd and 7de are the interfacial energies of the substrate-environment, 
droplet-substrate and droplet-environment interfaces, respectively. Therefore, 
this model can be used to simulate a case where the surface energies are 
spatially and/or temporally dependent on other fields, such as surface tem- 
perature or surface composition. This specific case applies when the wetted 
substrate behind the droplet have a higher interfacial energy than the pristine 
substrate, as in Ref. [80] . 



The evolution of the droplet surface is illustrated in Fig. [15} The droplet 
initially has the shape of a hemisphere, with a 90° contact angle with the sub- 
strate surface. The early evolution is marked by the evolution of the droplet 
shape as it relaxes to satisfy the contact-angle boundary condition, as seen 



in Fig. 15 a). The droplet then begins to accelerate. Once the contact angle 
reaches the prescribed value, it is maintained as the droplet moves toward 
the right; see Fig. [I5|(b). In the steady state, the droplet moves at constant 
speed without other effects present. Such motions of droplets have been ob- 
served and explained as a result of an unbalanced surface tension between the 
head portion (with a nonwetting surface) and tail portion (with a wetting sur- 
face) because of the resulting spatially varying composition and composition- 
dependent surface energy [80] . 



Figure 16 shows the relaxation of an initially hemispherical droplet on an 
irregular substrate surface in a 3D simulation. The contact-angle boundary 
condition imposed at the three-phase boundary is 135°. The computational 
box sizes are L^ = Ly = 120 and L^ = 80. As shown here, the droplet changes 
its shape to satisfy the imposed contact angle, and the droplet evolves into a 
shape for which the total surface energy is minimized. The behavior favoring 
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dewetting imposed by the contact angle [9 > 90°) is properly reflected in the 
lifting of the droplet, as shown in Figs. [I6|^a)-(c) and (d)-(f). During this 
relaxation process, the three-phase boundary moves toward the center as the 
droplet-substrate contacting area decreases, as shown in Fig. [l6|^a)-(c). This 
model and numerical method has been applied to simulate a nickel particle 
coarsening process in the complex channel within supporting porous ceramic 
microstructure (consisting of yttria-stabilized zirconia) in SOFC anodes, and 
to estimate the degradation of the anode material during SOFC operation 



6 Discussion and Conclusions 



In this work, we demonstrated a generalized formulation of the smoothed 
boundary method. This method allows Neumann, Dirichlet, or mixed bound- 
ary conditions to be imposed on a diffuse interface to solve partial differential 
equations within the region where the domain parameter ip uniformly equals 
1. The derivation of the method, as well as its implementation, is straight- 
forward. The method can be used to solve differential equations numerically 
without complicated and time-consuming structural meshing of the domain of 
interest, as the domain boundary is specifled by a spatially varying function. 
Instead, any grid system, including a regular Cartesian grid system, can be 
used with this method. 

This smoothed boundary approach is flexible in coupling multiple differen- 



tial equations. In Section |3.3[ we demonstrated how this method can be used 
to couple bulk diffusion with surface react ion- diffusion into a single equa- 
tion while the two equations serve as complementary boundary conditions. In 
principle, this method can be used to couple multiple differential equations in 
different regions deflned by different domain parameters. For example, if the 
physics within a domain deflned hj ipi = 1 are governed by a differential equa- 
tion Hi, the overall phenomenon will be then represented hj H = J^i'^tHi, 
where the subscript 'i' denotes the ith domain and Si V'* = 1 represents the en- 
tire computational box. When sharing the diffuse interfaces between domains, 
the physical quantities can be interconnected as boundary conditions for each 
equation in each domain. Therefore, this method could be used to simulate 
coupled multiphysical and/or multiple-domain problems such as fluid-solid 
interaction phenomena or diffusion in multi-material polycrystalline solids. 

We further demonstrated the capability of applying the smoothed boundary 



method to moving boundary problems in Section 5.2[ When the locations of 
domain boundaries are updated by a phase-fleld-type dynamics such that the 
domain parameters remain uniformly at 1 and on either side of the interface, 
the smoothed boundary method can be conveniently employed to solve partial 
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differential equations with moving boundaries. In addition to the phase-field- 
type dynamics, the smoothed boundary method is also applicable to moving 
boundary problems implementing the level set method P2|33f82] . with the 
domain parameter obtained simply by taking the hyperbolic tangent of the 
distance function. 

In addition to Neumann and Dirichlet boundary conditions, we also showed 
the capability of the smoothed boundary method for specifying contact angles 



between phase boundaries and domain boundaries (Sections 4.3, 5.4 and 5.5). 
This type of boundary condition is difficult to impose using conventional sharp 
interface models. 

Although the smoothed boundary method has many advantages, as shown 
in the results in Section |4] and in the derivations in Appendix [A], the nature 
of the diffuse interface inevitably introduces an error proportional to the in- 
terfacial thickness because we expand an originally zero-thickness boundary 
into a finite thickness interface. This spread of interface also leads to another 
error source depending on the resolution of the rapid transition of the do- 
main parameter across the interfacial region. When numerically solving the 
smoothed boundary formulated equations, properly capturing the gradient of 
the domain parameter across the interface becomes very important. Based on 
our experience, 3-5 grid points are necessary to properly resolve the diffuse 
interfaces while ensuring that the errors are well-controlled. Moreover, when 
solving time-dependent equations, one singularity occurs because of the terms 
l/?/' and l/ip"^ used to impose the Neumann and Dirichlet boundary condi- 
tions, respectively. In practice, a small value is necessarily added to ip to avoid 
singularity resulting from division by zero. In our simulations, the errors were 
quickly saturated when the value added to ip was selected to be smaller than 
1 X 10~^ when 3-5 grid spacings were used for the interfacial regions, which 
suggests it is unnecessary to select a smaller value for the singularity-control 
term. However, when solving time-independent equations, such as the mechan- 
ical equilibrium equation and the steady-state diffusion equation, there are no 
singular terms in the equations. The small additional term is then merely used 
to condition the matrix solver. In this case, it can be on the order of numerical 
precision, such as 1 x 10^^®. 

Based on the general nature of the derivation, the smoothed boundary method 
is applicable to generalized boundary conditions, including time-dependent 
boundary values important for simulating the evolution of many physical 
systems. Because the domain boundaries are not specifically defined in the 
smoothed boundary method, this method can be applied to almost any ge- 
ometry as long as it can be defined by the domain parameter. The developed 
method is thus a very powerful and convenient technique for solving differential 
equations in complex geometries that are often difficult and time-consuming to 
structurally mesh. As three-dimensional image-based calculations are increas- 
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ingly prevalent in scientific and engineering research fields J83f84 



where 

voxelated data from serial scanning or sectioning are often utilized and are dif- 
ficult to render as meshes, the smoothed boundary method is expected to be 
widely employed to simulate and study physics in complex geometries defined 
by 2D pixelated and 3D voxelated data with a simple process of smoothing 
the domain boundaries. 
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A Proof of Convergence for Neumann and Dirichlet Boundary 
Conditions for the Diffusion Equation 



To demonstrate that the smoothed boundary formulated diffusion equation 
satisfies the assigned Neumann boundary condition (specifying the boundary 
flux or normal gradient), we use the one-dimensional version of Eq. ^ without 
loss of generality for cases with higher dimensions. By reorganizing terms and 
integrating over the interfacial region, we obtain: 



«»+«/2 /dC 



ai-^/2 



^ 



\dt 



S ] dx = ijjD 



dC 
dx 



aiH/2 



a, -5/2 



ai+5/2 
a,-C/2 



dx 



DBjsidx 



(A.l) 



where a^ — ^/2 < x < Oj + ^/2 is the region of the interface and ^ is the 
thickness of the interface. Following Refs. [2|3ll20f23] . we introduce the mean 
value theorem of integrals, which states that for a continuous function g{x) 
there exists a constant value, /iq, such that viim.g{x) < Jp g{x)dx/{q — p) = 
Hq < Bi8ixg{x), where p < x < q. By eliminating the second terms on the 
right-hand sides of Eqs. (|6| and (A.l), the no-fiux boundary condition can be 
imposed {B^ = 0); the resulting equation is similar to those proposed in Refs. 
[T|2|3|20]|23] . However, here we retain the term to maintain the generality of 
the method. Therefore, the analysis presented herein leads to an extension of 
the original method that greatly extends its applicability. 



Because the function on the left-hand side of Eq. (A.l) is continuous and 



finite within the interfacial region, we can relate its value to the interfacial 
thickness by hg^, according to the mean value theorem of integrals. Using the 
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conditions that ip = 1 at x = Oj + ^/2 and ip = at x = Oj — ^/2, the first term 
on the right-hand side of Eq. ( |A.l ) is written as D{dC/dx)ai+^/2- Because 
\dip/dx\ = for X < Oj — ^/2 or x > a^ + ^/2, the bounds of the integral can 



be extended to — cxd and cxd. Therefore, we can rewrite Eq. (A.l) as 



ho^ = D 



dC 

dx 



ai+^/2 



dip 
dx 



DBjvdx 



(A.2) 



and by taking the hmit of this expression as ^ — )■ 0, we obtain: 



D 



dC 
dx 



/+00 
6{x — ai)DBi^dx = DBpf 



(A.3) 



where dC/dx\ai+^/2 — dC/dxla^ and hm^_i.o \dip/dx\ = 6{x — ai) when ifj takes 
the form of a hyperbohc tangent function and 5{x — ai) is the Dirac delta func- 
tion. The Dirac delta function has the property that j^^ 5{x — ai)f{x)dx = 



f{ai), providing the second equality in Eq. (A.3). Therefore, Eq. (A.3) clearly 



shows that the smoothed boundary method recovers the Neumann boundary 
condition at the boundary when the thickness of the diffuse boundary ap- 
proaches zero. This convergence has been observed for both stationary and 
moving boundaries [20|32|I33] . 

To demonstrate the convergence of the solution at the boundaries to the spec- 
ified boundary value, we again use a one-dimensional version of the smoothed 
boundary formulated equation. Integrating Eq. ([T]) over the interfacial region 
and reorganizing terms, we obtain: 



-«/2 



i^' 



,dC 

In 



^ 



d_ 

dx 



, dx 1 




dipC 
dx 



iP^S dx = - I ' D[^] ^-B 
Jai-i/2 \dx I dx 

(A.4) 

Similar to the derivation of Eq. (A.3), the left-hand side of Eq. (A.4) is propor- 



D 



dip 
dx 



dx. 



tional to the interfacial thickness and approaches zero in the limit of ^ — )■ 0. 
On the right-hand side of Eq. (A.4), the gradient of ip approaches the Dirac 
delta function, 6{x — ai), as the interface thickness approaches zero. Therefore, 
we can reduce Eq. (A.4) to limg_i.o /io^ = —D[d{ipC)/dx — B^dip/dx] in the 



limit .^ — 7- 0. By integrating over the interfacial region again, we obtain: 



hm ^ 
C^o D 



C 



ai+?/2 



«i+«/2 dip 

BD^::-dx, 

a,-i/2 dx 



(A.5) 



which gives C|a- = -BdU, in the limit of ^ — )■ because C\a, 



\im.^^Q{dip / dx) = 6{x 



+?/2 



C\a, and 



Therefore, the smoothed boundary formulation 



recovers the specified Dirichlet boundary condition: C = B^ at x 
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B Surface Laplacian Operator and Alternating Direction Line Re- 
laxation (ADLR) Method for Solving Coupled Surface Diffusion, 
Reaction and Bulk Diffusion Equation 



The surface gradient operator is defined by: 



' 1 — 77,1^1 — 721^2 —riin^ 



d - n(g)n)V 



V 



-77,2^1 1 - 77,2^2 -"-2^3 
-773771 -773772 1 - 773773 



\ 


d/dxi 




d/dx2 


/ 


d/dxs 



(B.l) 



where rii is the ith component of the inward unit normal vector (here, i = I, 
2 and 3, corresponding to the x, y and z directions, respectively). In tensor 
notation, this operator is written as Vs = rriijd/dxj. The repeated indices 
indicate summation over the index. The coefficients rriij are related to the 
surface unit normal by mn = 1 — 77i77i, 77122 = 1 — ''^2^2, ""^33 = 1 ~ ''^3^3; 
77^12 = 77^21 = -nin2, 777-13 = "^-si = -nim, and 77^23 = ^■32 = -^2't-3- The 
surface Laplacian operator is defined by the surface divergence of the surface 
gradient: 

The scalar surface Laplacian is a sum of nine second-order-partial-differential- 
operator terms (where j = k) and 18 mixed-partial (cross) terms of the dif- 
ferential operator (where j ^ k). 



To solve Eq. (26), which includes the surface Laplacian operator, we use an 



ADLR method in a second-order central difference scheme. We first separate 



the surface Laplacian operator in Eq. (B.2) into a term involving pure second- 
order partial derivatives ("diagonal") and another containing mixed partials 



cross 



V^ = V^- + V^ , where: 
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(B.3) 
Thus, Eq. (|26|) is rearranged to [V ■ (^I^bV) - |V^|(fi; - /^sVLg)]C' = {iuip - 
|V'?/'|-DsV^ross)C', where the three axial components of the V- (V'-DbV) operator 
and the "diagonal" terms in the V^jag operator can be discretized by central 
difference schemes similar to: 



dx\ dx I Ax 



Wi+l/2,j,k T Wi-l/2j,k — 



Ax 



Ax 



(B.4) 
and ipi+i/2,j,k = ('^i+i,j,fc+V^j,j,fc)/2. The discretization scheme along the y and z 
axes can be similarly obtained. Therefore, the operator [V- (■?/'-DbV) — | V'?/'| (ft — 
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DgV^iag)] can be lumped into an equivalent Helmholtz operator containing 
only second-order partial derivatives of the neighboring six points along the 
X, y and z axes and a coefficient term IV-j/'Ik at the center point. Additionally, 
the 18 terms in the V^oss operator can be calculated using a discretization 
scheme. For example, the d{mdC/dy)dx can be discretized as: 



dx\ dy 



a 



2Aa; 



m 



i+lj,k' 



i+l,j+l,k 



2Ay 



(^i+l,j-l,k (^i-l,j+l,k 



a 



i—l,j—l,k 



2Ay 



(B.5) 

and similarly for other components. As a result, the equation is represented 
by LC = §, where £ is a linear Helmholtz-like operator and S is calculated 
as {iuip — I VV'l-DsV^j.QgJC*. This equation can be solved using an ADLR solver 
[86 .87] by decomposing the Helmholtz operator into the three axial directions: 



P ^(n+l/3) ,,-, ^(n+1/3) , p ^(n+1/3) _ p{n) n ^(n) n ^{n) 

'^i+l,j,k^i+lj^k ~ ^i,j,k^i,j,k + '^i-^,j,k^i-l,j,k - ^ " ^yy^yy ~ ^zz'^zz i 

(B.6a) 



r /^(n+2/3) r.,-. ^yiL-r^/o/.p 

'^ij+l,k'^ij+l,k ~ ^^i,j,k'-^i,j,k +'^ 



^(n+2/3) ^ ^{n+2/3) _ o(„) n /=y(n+l/3) r /=Y(n+l/3) 



■'xx^xx 



L 



i,j,k- 



-iC, 



(n+1) 



i,j,fc+l ^i,j,k^i,j,k 



(n+1) ^ ^(n+1) 



^i,j,k- 



''i,j,k—l 



zz^zz 

(B.6b) 

«(")-£ Ci^+V3)_r -^(n+2/3) 
" ^xx^xx ^yy^yy ' 

(B.6c) 

^i,j+l,k(-^i,j+l,k + 

'^i,j-i,kCi,j-i,k, '^zzCzz = ^i,j,k+iCi^j^k+i + ^i,j,k-iCi,j,k-i, and the supcrscript 
(n) denotes the nth iterative step. Within each iterative step, we ffist solve 



where •Lxx^-'xx — -^i+lj.fct-^i+lj.fc + ^i-l j,k^i-l,j,ki yyy^yy 



along the x direction using Eq. (B.6a), for which a simple tridiagonal matrix 



solver is employed for each column. Similarly, Eqs. (B.6b) and (B.6c) are solved 



along the y and z directions, respectively, with the updated value on the right- 
hand sides. The above procedure is repeated until the solution converges to 
its equilibrium value. 



For the simulations of oxygen-vacancy diffusion in a cylinder in Section 4.2 



we consider a cylindrical symmetry for the differential operator, in which the 



dimensions reduce to effectively 2D, such that the "bulk" term in Eq. (12) 
becomes: 



r or \ or I 



with only components in the radial and axial directions. Here, we have set Db 
at 1 for clarify of the derivation. The ffist term on the right-hand side can be 
rewritten and discretized using the central difference scheme as: 



d\ rip — 

d{r'^) \ dr 



'i+l/2j ' i-l/2,j 



^i+l/2,j'^i+l/2,i 



c,. 



i+l,i 



^ij 



^i+l,j ^i,j 



- ri-i/2,j'ipi-i/2,j 
(B.8) 



(-^i.-j (-^', 



'^J 



i-i,j 



^i,j ^i-l,j 



where the subscript i and j denote the ith and jth grid points in the radial 
and axial directions, respectively. If the radial grid spacing is selected to be 
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uniform, r-j+i ,• - r.^ = ru - ri^^ 



Ar. The second term on the right-hand 



side of Eq. (B.7) can be discretized using the scheme provided in Eq. (B.4). 



The surface Laplacian with cyhndrical symmetry is given as: 



„2^ 1 d ( dC\ d ( dC\ d ( dC\ d ( dC\ 

V^C =mrr-^\ rnirr^- +m„,— m^^-- +mr2-- m„-- +^^^--771^^-- + 
r or \ or J or \ oz J oz\ or J oz\ oz J 



m. 



l_d_ 
r dr 



rm. 



dr I 



+ m. 



d 



dr \ ^^ dz ) 



+ m. 



d 



m 



dC\ 



d 



or J 

(B.9) 



dC 



; O \ ■ '"Zr n I I ' ■ "ZZ n I " "2:2 n I 5 

OZ \ or J OZ \ OZ 



where rrir 



X I Lv' I Lv' t, llli'py 



m. 



-UrUz, and ruz 



UzUz- The 



"diagonal" terms along the radial and axial directions can be discretized in a 



manner similar to Eqs. (B.8) and (B.4), respectively, and the "cross" terms 



can be discretized as in Eq. (B.5 ) for solving Eq. ( 12 ) with an additional factor 
of m. 



C Derivation of the Mechanical Equilibrium Equation 



To perform the smoothed boundary formulation on the tensorial mechanical 



equilibrium equation, we multiply Eq. (14) by ip and use the mathematical 



identity ijj{dHij/dxj) = d{ipHij)/dxj — {dip /dxj)Hij to obtain: 



_d_ 

dxi 



^a 



ijkl 



1 f duk du 



2 V dxi dx 



dtp Ifduk dui^ 
dxj ^'^ 2 1 dxi dxk 



- — I ippCijki^ki I — ^ — pCijkiSki- 
(C.l) 



By collecting the terms associated with dijj/dxj, we obtain Eq. (15) as given 
in Section Em 



To impose a Dirichlet boundary condition (a specified displacement) on the 



mechanical equilibrium equation, we multiply the left-hand side of Eq. (14) 
by ip"^ to obtain: 



2 <9 „ 1 fduk , dui 

ip - — Ci- 



dxi 



^ijkl 



2 V dxi dx 



ij 



_d_ 
dx-i 



'^''^2\dxi dxk 



dip I fduk , dui 



dx^ 



2 y dxi dxk 
(C.2) 



where the second term on the right-hand side can be replaced by: 



dxj ""^ 2\ dxi dxk j dxj ^^ 2 



d{ipUk) d{%pui) 



dxi 



dxi 



dip I ( dip dip^ 

■^—CijM- \Uk^ h ui- — 

dXj 2 \^ dxi dxk 
(C.3) 
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according to the product rule: 



dil) 
dxi 



a 



ijkl 



d{ipUk) d{i!Ui) 



dxi 



dxu 



dip 

dXn 



a 



ijkl 



duk , d^ 
. dxi dxi ^ 



( , dui dih 



dxi 



' dxi 



dip fduk , dui 
--W^^<-"ijki\ - — \- 



dxj 



dxi 



Thereby, we obtain Eq. (18) in Section 3.4 



d^) f dip dip' 

dxkj dxj \ dxi dxk^ 

(C.4) 



D Relation Used in the Derivation of Contact Angle Boundary 
Condition 



Here, we multiply the equilibrium criterion of the phase field model by V0 to 
obtain: 



df 



V0 - (e2v20)V0 



V/ - -V{V<pr 



0. 



(D.l) 



Integrating the above, we obtain / — e^|V0p/2 = Ci, where Ci is a constant 
of integration. In the phase field model, the order parameter remains at a 
uniform value in the bulk away from the interface; thus giving |V0| = in 
the bulk. Therefore, ci is equal to the bulk value of /. For convenience, we 
have taken the free energy at the bulk values to be zero, and therefore ci = 0, 
leading to V(p = V27/e. However, the choice of the free energy value at the 
bulk is arbitrary, and therefore does not affect the result of the calculation as 



long as it is taken into account by replacing / appearing in Eq. (20 ) by / — Ci 



E Smoothing Voxelated Data Using a Distance Function 



The experimentally obtained microstructure is typically provided in a form 
of a 3D array containing voxels of different values indicating different phases. 
To incorporate the voxelated data into the smoothed boundary formulation, 
we must convert the discrete voxelated array into a domain parameter profile 
that continuously transitions from one phase to another. Here, we employ 
the distance function method commonly used for initialization in the level 
set method |12|55] . First, we construct the sign function by assigning positive 



and negative values to the voxels in the solid and pore phases, respectively: 
Sgni^s.) = 1 for the solid phase and Sgn{'x) = — 1 for the pore phase, where 
X is the position of a voxel. The distance function indicating the distance 
between the center of a voxel and the sohd-pore interface is calculated by 
evolving the time-dependent equation dip{'x,t)/dt = Sgn{x){l — |Vv9(x, t)|) 
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to its equilibrium. This process is numerically implemented by: 

r(^(x,t) + At[5(7r2(x)(l-|V<^(x,t)|)] if (^(x,t + At)-(^(x,t)>0 
V?(x, t + At) = < 

I !f{x,t) + At[Sgn{x)v] if ¥?(x, t + At) ■ (/?(x, t) < 

(E.l) 
where t; is a small nonzero value. The second case above prevents interfaces 
from moving more than one grid spacing by requiring that the sign of the 
function remain the same as the initial value. The absolute value of the gra- 
dient of the distance function is calculated using a Godunov upwind scheme 



\Vipij^k\ = [max(max(L)+v9jj-fc,0)^max(-Z:'^, (/}ij^fc,0)^) + 

max(max(£)+v9ij- fc, 0)^, max{-D~ipij^k, 0)^)+ (E.2) 

max(max(D+V9ij- fc, 0)^ max{-D^ipij^k, Of)Y^^, 

where i, j and k are the indices of the grid points along the x, y and z axes, 
respectively, and: 

Dy(pi,j,k = {^i,j,k - (pi,j-i,k)/Ay, Dyipij^k = {^i,j+i,k - (pi,j,k)/Ay, (E.3) 

Dt'^i,j,k = {'^i,j,k - '^i,j,k-i)/^z, D'ipij^k = (v'ij.fc+i - Vij,k)/^z. 

In practice, for the smoothed boundary method, the function must take the 
form of the distance function only near the interfacial regions, and therefore 
the convergence condition can be placed in these regions only (and not in the 
bulk far from interfaces) as long as the values in the bulk are sufficiently large 
in magnitude. 

From the distance function, we obtain a domain parameter based on the ex- 
perimentally acquired voxelated data by taking the hyperbolic tangent of the 
distance function, ^(x) = {1 + tanh[99(x)/^]}/2, where ip = 0.5 coincides the 
location of the zero level set [ip = 0), ip = 1 in the solid, ip = in the pore, 
and the value of ( controls the thickness of the interface. 



F Projection Method 



To simulate Kirkendall-effect-induced deformation, we model the solid diffu- 
sion couple as a very viscous fluid that deforms in a quasi-steady-state manner, 
namely, creep flow. In contrast, the environmental phase surrounding the solid 
is treated as a nearly inviscid fluid. A simple way to implement this model 
is to deflne the viscosity coefficient as ri{ip) = fjip + v, where f] is a constant 
viscosity coefficient for the solid phase and i; <^ r^ is a small value used to 
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avoid numerical instability. To solve the velocity field with a variable viscosity 
coefficient, we adopt a projection method jHSES], in which the divergence of 



the viscous stress tensor is decomposed into a linear part and a residual part, 



givmg: 



V ■ r/[Vv + (Vv)^ ] = AV ■ [Vv + (Vv) 



T^ 



(F.l) 



where A is a constant scalar numerical parameter for the scheme (normally 
set between 0.5?7 and ff) and r^ is a vector residual. Using the identity that 
V-[Vv + (Vv)^] = VV + V(V-v), where V-v = gy, and V^v = d^Vi/dxjdxj 
is a vector containing the Laplacian of each velocity component, we rewrite 



Eq. (28a) as: 



VP + AVV + rv + V A 



|^y + ^/iVV' = 0. (F.2) 



By taking the divergence of Eq. (F.2), applying the relation V ■ (V^v) = 
V^(V ■ v) = V'^Qv and rearranging the terms, we obtain the Poisson equation 
of the scalar pressure field, which serves as one of the two equations for the 
iterative scheme: 



y2p(n) _ y . ^^(n-l) ^ y2 2A 



|^y + ^V-(/iV^) 



(F.3) 



where the superscript (n) denotes the nth iterative step. The second and third 
terms on the right-hand side are fixed during an evolution time step while the 
values of the velocity and pressure are updated during iteration. 



For the velocity field, we reorganize Eq. (F.2) to obtain the Poisson equation 



for the velocity component in each coordinate direction: 



V2v(") 



A 



Vp(") 



("-!)_ V A 



9v 



a 



iiVij 



0. (F.4) 



The residual vector Vy, is calculated using Eq. (F.l) and is updated during the 



(n) 



iteration: Tv^"^ = V ■ r^fVv^") + (Vv^"))^] - AV^v^") - AVgv- The Poisson 
equations can be solved using an ADLR method similar to that described 
in Appendix |B} except that the Helmholtz operator is replaced here by a 
Laplacian operator. Within each time step for the deformation (Eq. (28c)) 



and diffusion (Eq. (30)) of the diffusion couple, the pressure and velocity fields 



are solved iteratively until the values of pressure and the velocity components 
converges. The convergence criteria is set to be e < 1 x 10^^, where e is the 
relative error taken by dividing the root-mean-square difference between the 
values of two consecutive iterative steps by the average magnitude of the values 
in the previous step. The velocity field is then substituted into the advective 
terms in the order parameter and concentration evolution equations. 
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G Implicit Time Scheme 



Equation (32) contains the coupled diffusion equations for two species and 
is constrained by a very small time step for the Euler explicit time scheme 
because of the large diffusivity, Dyv (nearly 10^ times larger than D^^). 
Thus, here we use a semi-implicit time scheme, similar to that presented in 
Ref. [72j, to solve the vacancy diffusion equation and to significantly enhance 
numerical efficiency. In the time-discretized form, the scheme is given by: 

At V^(") ^ 



(G.l) 

where the superscript (n) denotes the nth time step and x is a weighting 
factor that can be optimized to increase numerical stability. The diffusivities, 
Dvv and Dvb, are mole-fraction dependent quantities. The average diffusivity, 
Dyy, is calculated from Dyy over the solid region in which diffusion occurs. 



The diffusion equation for B atoms, Eq. (32b ), and the Cahn-Hilliard equation 



Eq. (33), are solved using the Euler explicit time scheme, as the diffusivities 



and mobilities for these equations are much smaller. 



H Mechanical Equilibrium Equation Solver 



Here, we expand the generalized mechanical equilibrium equation, Eq. (17), 
for a linear elastic and isotropic solid. In this case, the components of the 
elastic constant tensor are expressed by: 

All = C*llll = C*2222 = C*33335 (H.la) 

(H.lb) 



Al2 


= 


C1122 — 


C22II — 


C2233 — 


C3322 — 


C33II = 


Cll33; 




A44 =Ci212 


= 


C122I 


= 


C2II2 


= 


C212I 


= C*2323 


= 


C2332 






— C*3223 


= 


C3232 


= 


C1313 


= 


C133I 


= Czuz 


= 


C313I 



(H.lc) 

The remaining elastic constant components vanish. For an isotropic solid, the 
Young's moduli are related to Lame constants hy E = Ai2(l + z^)(l — 2z/)/i/, 
where v is the Poisson's ratio; the shear modulus is given by A44 = Ai2(l — 
2z/)/2z/, and the elastic constant An is given by An = A12 + 2A44. The pair, 
A12 and z/, forms the set of Lame constants. 

We use coordinate notation to replace the indices i = 1, 2 and 3 with x, 
y and z, respectively. With a traction-free boundary condition on the solid 
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surface {Ni = 0), the smoothed boundary formulated mechanical equilibrium 



equation, Eq. (17), can be written out for the x, y and z directions as: 



d_ 
dx 
d_ 
dx 



(du\ 
^\dx) 


d 
dy 


1 


'du\ 

.9y) 


d 

+ dz 


ipX^ 


'duV 
^dz) 


d 

dx 


(dv 
\dy 


dw\ 
dz)\ 


d 

dy 


i/jXu 


(dvV 
[dx)_ 


d 

dz 


^A44 


(dwV 
[dx)_ 



[V^p(Aii + 2Ai2)]- 



(H.2a) 



d_ 
dx 
d_ 
dx 




+ 



d_ 

dy 

d_ 

dy 




[#(Aii + 2Ai2)]- 



(H.2b) 



d_ 

dx 

d_ 

dx 




[MA11 + 2A12)]- 



(H.2c) 



where m, v and w are the displacements along the x, y and z axes, respectively. 
Here, we provide the long form of the mechanical equilibrium equation for the 



sake of clarity for the readers. To solve Eq. (H.2), we use an ADLR method 



similar to that provided in Appendix |Bj We keep the "diagonal" terms, such 
as d{'ipXd/dx)/dx, d{il)Xd/dy)/dy, and d{il)Xd/dz)/dz on the left-hand sides, 
and move the "cross" terms of the differential operator to the right-hand sides. 
As a result, we obtain three equations with second-order-partial-differential 
operators for the three displacement components. The "diagonal" terms on 



the left-hand sides can be discretized using the scheme in Eq. (B.4). The 
"cross" terms moved to the right-hand sides can be calculated in a scheme 



similar to Eq. (B.5). The ADLR method shown in Eq. (B.6) is employed to 



solve the displacements, and the solutions are updated and iterated until all 
the displacement components reach their equilibrium values. 



For the case in Section 5.3, the solid phase {ip = ipi+ V'2) includes two differ- 
ent materials: GDC {ipi) and LSC {'ipi)- Therefore, we smoothly interpolate 
material properties appearing in the differential operators: ifjXij = ipiX^^ 
ipiXif" . Similarly, the body force term is interpolated by ipp{Xii + 2A 

jiDC(\GDC I o\GDC\ I ,/, „LSC(\LSC 



GDC 



12 



7/^iP^^^'(Afi^^ + 2Af2^^) +V^2P^^^(Aff^' + 2Af2^'^'), where p = a AT is the ther- 
mal expansion. The thermal stress is calculated according to Hooke's Law, 
which is written with the domain-parameter-interpolated elastic constants and 
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thermal expansions as: 

(^ij = {Wl<^ijkl +^2Cij-fc, )--— + -— -(^ip Cijkl +W2P Cijkl >Okl- 



2 V dxi dxk , 



(H.3) 
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Case 


1 


2 


3 


4 




Ax = 2.5 X 10-2 
V = lx 10-^ 


C = 1.145AX 
v = lx 10-^ 


C = 5.73 X 10-2 
V = lx 10-^ 


Ax = 2.5 X 10-2 
C = 2.86 X 10-2 


c 


e 


Ax 


e 


Ax 


e 


V 


e 


a 


1.43 X 10-2 


2.74 X 10-1 


1.25 X 10-2 


3.93 X 10-1 


1.25 X 10-2 


1.75 X 10-3 


1.0 X 10-2 


7.75 X 10-3 


b 


2.86 X 10-2 


*7.88 X 10-1 


2.50 X 10-2 


*7.88 X lO-i 


2.50 X 10-2 


«1.72x 10-3 


1.0 X 10-3 


1.39 X 10-3 


c 


5.73 X 10-2 


*1.72 X 10-3 


5.00 X 10-2 


■^1.58 X 10-3 


5.00 X 10-2 


^1.58 X 10-3 


1.0 X IQ-' 


7.93 X lO-i 


d 


1.15 X 10-1 


3.53 X 10-'^ 


1.00 X 10-1 


3.20 X 10-3 


1.00 X 10-1 


1.16 X 10-3 


1.0 X 10-^ 


*7.88 X 10-1 


e 


2.29 X 10-1 


7.20 X 10-'^ 


2.00 X 10-1 


6.54 X 10-3 


2.00 X 10-1 


7.53 X 10-'' 


1.0 X 10-'-' 


7.88 X lO-i 


f 


4.58 X 10-1 


1.49 X 10-2 


4.00 X 10-1 


1.39 X 10-2 


4.00 X 10-1 


unstable 


1.0 X 10-11 


7.88 X lO-i 



Table 1 

Relative errors, e, for the ID smoothed boundary diffusion equations to the analyt- 
ical solutions with various parameters. Each of the markers, *, o and <i, denotes the 
identical result from a set of parameters. 





K = 2.1 


K = 20 


«; = 50 


K= 100 


thin-interface 
Co = 0.075 


e 


7.99 X 10-1 


e 


2.26 X 10-3 


e 


2.46 X 10-3 


e 


7.26 X 10-3 


eb 


7.99 X 10-1 


eb 


2.23 X 10-3 


eb 


2.39 X 10-3 


eb 


7.37 X 10-3 


es 


8.03 X 10-1 


es 


3.02 X 10-3 


es 


4.02 X 10-3 


es 


2.63 X 10-3 


medium- inter face 
Co = 0.149 


e 


1.08 X 10-3 


e 


3.06 X 10-3 


e 


8.32 X 10-3 


e 


2.74 X 10-2 


eb 


1.04 X 10-3 


eb 


2.89 X 10-3 


eb 


8.51 X 10-3 


eb 


2.84 X 10-2 


es 


1.50 X 10-3 


es 


4.83 X 10-3 


es 


5.12 X 10-3 


es 


1.86 X 10-3 


thick- inter face 
Co = 0.292 


e 


1.81 X 10-3 


e 


1.08 X 10-2 


e 


2.90 X 10-2 


e 


7.34 X 10-2 


eb 


1.63 X 10-3 


eb 


1.13 X 10-2 


eb 


3.11 X 10-2 


eb 


7.89 X 10-2 


es 


2.69 X 10-3 


es 


6.78 X 10-3 


es 


5.75 X 10-3 


es 


1.06 X 10-3 



Table 2 

Relative errors for the coupled surface reaction-diffusion and bulk diffusion model 

in a cylinder, where e denotes the overall error, Cb denotes the bulk error excluding 

the surface points, and Cg denotes the surface error calculated only with the surface 

points. 
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Contact angles for Allen-Calm equation 


c 


5^ = 1.0607 


5^ = 1.4142 


5^ = 1.7678 


S^ = 2.1213 


5^ = 2.8284 


0.75 


0.5050 (59.67°) 


0.5048 (59.68°) 


0.4965 (60.23°) 


0.5001 (59.99°) 


0.5004 (59.97°) 


1.00 


0.4900 (60.66°) 


0.5039 (59.74°) 


0.4966 (60.22°) 


0.4982 (60.12°) 


0.4956 (60.29°) 


1.50 


0.4865 (60.89°) 


0.4962 (60.25°) 


0.4927 (60.48°) 


0.4938 (60.41°) 


0.4927 (60.48°) 


2.00 


0.4886 (60.75°) 


0.4918 (60.54°) 


0.4883 (60.77°) 


0.4901 (60.65°) 


0.4901 (60.65°) 


4.00 


0.4825 (61.15°) 


0.4782 (61.43°) 


0.4783 (61.43°) 


0.4795 (61.35°) 


0.4790 (61.38°) 






Contact angles 


for Cahn-Hilliard « 


quation 




C 


S^ = 1.0607 


S^ = 1.4142 


S^ = 1.7678 


5^ = 2.1213 


S^ = 2.8284 


0.75 


-0.4831 (118.89°) 


-0.5003 (120.02°) 


-0.4937 (119.59°) 


-0.4979 (119.86°) 


-0.4931 (119.54°) 


1.00 


-0.4923 (119.49°) 


-0.4926 (119.51°) 


-0.4897 (119.32°) 


-0.4929 (119.53°) 


-0.4881 (119.21°) 


1.50 


-0.4841 (118.95°) 


-0.4871 (119.15°) 


-0.4868 (119.13°) 


-0.4890 (119.28°) 


-0.4867 (119.13°) 


2.00 


-0.4639 (117.64°) 


-0.4857 (119.06°) 


-0.4861 (119.09°) 


-0.4862 (119.09°) 


-0.4853 (119.13°) 


4.00 


-0.4372 (115.93°) 


-0.4713 (118.12°) 


-0.4752 (118.37°) 


-0.4745 (118.33°) 


-0.4752 (118.37°) 




Dr 


ler parameter conservation for Cahn-Hilliard equation 




C 


S^ = 1.0607 


S^ = 1.4142 


5^ = 1.7678 


S^ = 2.1213 


S^ = 2.8284 


0.75 


0.9929 


0.9972 


0.9979 


0.9982 


0.9986 


1.00 


0.9930 


0.9973 


0.9979 


0.9982 


0.9986 


1.50 


0.9933 


0.9974 


0.9980 


0.9983 


0.9987 


2.00 


0.9976 


0.9991 


0.9993 


0.9994 


0.9996 


4.00 


0.9982 


0.9993 


0.9995 


0.9996 


0.9997 



Table 3 

Contact angle results from the validation simulations using the Allen-Cahn equation 
and the Cahn-Hilliard equation. In addition, a measure of the order parameter con- 
servation, evaluated by J 'il}(p{tss)dQ/ J ip(j){t = 0)dO,, where tgs is the time required 
to reach steady-state conditions and Q is the computational domain, are presented 
for simulations with the Cahn-Hilliard equation. 




Fig. 1. (a) Conventional sharp interface description of a domain bound by a ze- 
ro-thickness boundary, (b) Diffuse interface domain and boundary defined by a 
continuous domain parameter, ip. (c) Inward normal vectors defined by 'Vtp plotted 
for the square region in (b). 
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Fig. 2. Demonstration of the smoothed boundary method on the ID diffusion equa- 
tion. The red hne with circular markers is the domain parameter, and the blue 
lines are the concentration profiles at different times. The Neumann and Dirichlet 
boundary conditions are imposed at the right and left boundaries, respectively. The 
black lines are the corresponding analytical solutions of the sharp interface version 
of the diffusion equation. 
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Fig. 3. (a) Equilibrium concentrations in domains of two different interfacial thick- 
nesses, corresponding to Case If {( = 4.58 x 10^^) and Case lb ((" = 2.86 x 10^^) in 
Tablefl] The black line is the analytical solution, (b) Relative errors for the smoothed 
boundary solutions during diffusion for different (" values. Magnified views of (a) at 
the (c) left, and (d) right interfaces. 
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(a) |-e-V^2/ C'2f^iJ2b -C2b-Car,a\ (b) \-e2a-e2b e2c-e2d e2e-^2f\ 




250 500 750 1000 



Fig. 4. (a) Equilibrium concentrations for Cases 2f and 2b in Table [!} The black 
line is the analytical solution, (b) Relative errors for the smoothed boundary so- 
lutions during diffusion for the parameters in Case 2 in Table [ll (c) Equilibrium 
concentrations for Cases 3e and 3a. (d) Relative errors for Case 3 in Table [I] 
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Fig. 5. Steady-state concentration profiles for D]^ = 1, Dg = 10. The left column 
is for K = 2.1, and the right column is for k = 50. The top row shows the sharp 
interface solution; the middle row is the thin-interface smoothed boundary method; 
and the bottom row is the thick-interface smoothed boundary method results with 
interfacial thickness four times larger than those in the middle row. The top regions 
of constant blue color in (c)-(f ) represent the areas outside of the solid, whereas the 
solid white lines indicate the solid cylinder surface. 




4 xlO Q 0.02 0.04 0.06 




Fig. 6. Profiles of relative error between the sharp interface results and the smoothed 
boundary results having thin, medium and thick interfaces from the top to the 
bottom: (a), (c) and (e) have D^, = 1, Dg = 10, and n = 2.1, and are compared to 
the sharp interface result shown in Fig. [Sjl^a); (b), (d) and (f) have Db = 1, -Dg = 10, 
and K = 50, and are compared to the sharp interface result in Fig. ^ 
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Fig. 7. (a) Contour lines of ^ = 0.5 at various times for Allen-Cahn-type evolu- 
tion equation with a 60° contact-angle boundary condition, (b) Magnified view of 
the order parameter and (c) cos 9 profiles near the three-phase boundary (further 
magnified), corresponding to t = 2.7 x 10^ in (a), (d) Contour lines of tp = 0.5 at 
various times for Cahn-Hilliard-type evolution equation with a 120° contact-angle 
boundary condition, (e) Magnified view of the order parameter and (e) cos 9 profiles 
at the three-phase boundary (further magnified), corresponding to t = 3.0 x 10^ in 
(d). The cosine values of the imposed contact angles are 0.5 and -0.5 for (a)-(c) and 
(d)-(f), respectively. The circular markers in (c) and (f) denote the grid points in 
the range of 0.1 < ip < 0.9 and 0.1 < (/> < 0.9. The order parameters in the re- 
gion of V' < 0.5 have no physical significance. For the Cahn-Hilliard case, the order 
parameter is conserved in the region of '0 > 0.5. 
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Fig. 8. Simulation results of dimensionless steady-state oxygen-vacancy concentra- 
tions under DC loading for D\) = 1 using an experimentally obtained SOFC cathode 
complex microstructure as the input geometry: (a) k = 0.1 and D^ = 0; (b) k = 2.1 
and Ds = 10. 




Fig. 9. Simulation results of dimensionless steady-state oxygen-vacancy concentra- 
tions under AC loading for Z?b = 1, k = 2.1, and Ds = 10 using an experimentally 
obtained SOFC cathode complex microstructure as the input geometry. In (a) and 
(d), the real and imaginary parts for cj = 1.5 are shown, respectively. Similarly, (b) 
and (e) display the real and imaginary parts for uj = 51.5, respectively. Shown in 
(c) and (f) are the magnified views of (b) and (e), respectively. 
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Fig. 10. (a) Solid phase containing cathode (yellow) and electrolyte (semitransparent 
cyan) in a solid oxide fuel cell material, (b) Interfaces defined by y^|V'(/'i||V'02| > 0.2 
(semitransparent cyan) and \/|V^02|TVV^ > 0.2 (semitransparent purple). Dimen- 
sionless steady-state oxygen-vacancy concentration with the boundary condition 
C = 1 at the electrolyte-cathode interface for Db = 1: (c) k = 0.1, and D^ = 0; (d) 
K = 2.1, and D^ = 10. 
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Fig. 11. Left column: Fast diffuser mole fraction profiles (normalized to the lattice 
site density) recorded at dimensionless time of t = 0, 3.96 x 10 , 5.09 x 10^, and 
4.32 X 10^. Right column: Velocity fields corresponding to the mole fraction profile 
on the left. Black and gray arrows denote the flows inside and outside the material, 
respectively. The flow outside of the material has no physical significance for the 
shape change. Simulation parameters were 77 = 1 x 10^, Ca = 1 x 10^, e = 1 and 
M = 1.25 X 10~^. The fast diffuser hop frequency is four times larger than that of 
the slow diffuser. 
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Fig. 12. Top row, (a)-(d): snapshots of the fast diffuser mole fraction recorded at 
four different dimensionless times, t = 0, 2.79 x 10^, 2.62 x 10^, and 3.74 x 10^, 
respectively. Bottom row, (e)-(h): snapshots of vacancy mole fraction normalized 
to the equilibrium value, corresponding to (a)-(d). The solid white contour lines 
indicate the locations of the rod and void surfaces. The fast diffuser hop frequency 
is four times larger than that of the slow diffuser. 




Fig. 13. The mean stresses resulting from thermal expansion in (a) the entire solid 
phase, (b) the cathode phase, and (c) the electrolyte phase after rotating the volume 
180° around the z-axis. The unit of stress is Pa. 
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Fig. 14. (a) The evolution of an evaporating droplet on a rough surface (dashed 
line) governed by Allen-Cahn dynamics. The contact angle between the droplet 
and the surface is imposed at 135°. Solid curves of various colors represent the 
profile of the droplet at different times. The outermost blue line represents the 
initial state, and the innermost red line represents the final state (recorded before 
complete evaporation in the simulation); the lines are plotted at time intervals of 
270 dimensionless time units. The velocity of the three-phase boundary is greatly 
affected by the surface profile, (b) A magnified view of the three-phase boundary, 
showing that the contact angle is accurately set. The angle made by the thin black 
lines is 135°. (c) The order and domain parameters are shown to illustrate the diffuse 
natures of the interface and boundary (recorded at t = 270). 
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Fig. 15. A self-propelled droplet driven by unbalanced interfacial energies. The evo- 
lution is modeled by the Cahn-Hilliard equation with two different contact-angle 
boundary conditions on each side of the droplet, (a) The droplet shape changes 
during the relaxation period. The color contours are plotted at time intervals of 
2 X 10^ in dimensionless time units, (b) Droplet motion along the substrate surface. 
The color contours are plotted at time intervals of 1 x 10"' in dimensionless time 
units. The droplet moves at constant speed at steady state. 




Fig. 16. A droplet relaxing toward its equilibrium shape. The evolution is modeled by 
the Cahn-Hilliard equation with a contact angle of 135° to the irregular substrate 
surface: (a) initial {t = 0), (b) intermediate (f = 3 x 10^), and (c) equilibrium 
(t = 2.35 X 10^) states. The three-phase boundaries are delineated in red. Side 
views of the droplet are shown in (d), (e) and (f), corresponding to (a), (b) and (c), 
respectively. 
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